Mohammad Hossein Shaterzadeh-Yazdi

Análise de contato entre dois corposelásticos usando

o Método dos Elementos de Contorno

94/2015

CAMPINAS2015

i

Dedico este trabalho aos meus queridos pais;Meus primeiros e eternos professores.

vii

Agradecimentos

A Deus por minha vida, família e amigos.

Aos meus queridos pais, Azam e Mohsen, e às minhas irmãs pelo amor, incentivo e apoio incondi-cional.

Ao meu orientador, professor Dr. Paulo Sollero, pela oportunidade e apoio e também pelos conse-lhos, acadêmicos e não acadêmicos, que vou levar para resto da minha vida.

Ao meu co-orientador, professor Dr. Éder Lima de Albuquerque, pela motivação e indicação pararealizar o meu mestrado e também pelo apoio e pelas discussões construtivas na construção doprojeto.

Aos colegas e companheiros do laboratório Sollero-Pavanello, em especial René, Andrés e Kevinpela amizade e apoio.

À comunidade iraniana de Campinas e ao pessoal da casa CEU pela convivência e por tornarpossível morar longe da família e sentir menos falta de casa. Em especial agradeço ao Samu pelaamizade e ajuda na correção deste texto.

Aos professores e funcionários da Faculdade de Engenharia Mecânica da Universidade Estadual deCampinas- UNICAMP pela infraestrutura fornecida, apoio acadêmico e simpatia.

Ao Conselho Nacional de Desenvolvimento Científico e Tecnológico - CNPq pelo suporte finan-ceiro e por tornar possível a realização deste projeto.

Aos amigos que, mesmo de longe, me apoiam e não me deixam sozinho e a todos os que direta ouindiretamente fizeram parte da minha formação.

ix

Procure ser um homem de valor, em vez deser um homem de sucesso.

Albert Einstein

xi

Resumo

SHATERZADEH-YAZDI, Mohammad Hossein. Análise de contato entre dois corposelásticos usando o Método dos Elementos de Contorno. 2015. 102p. Tese (Mestrado). Faculdadede Engenharia Mecânica, Universidade Estadual de Campinas, Campinas.

Em problemas de contato mecânico entre dois corpos elásticos, o cálculo de tensões e defor-mações dos componentes é de grande importância. Em casos particulares os corpos estão sujeitosa cargas normal e tangencial na presença de atrito, o qual aumenta a complexidade do problema.O estudo do fenômeno e a modelagem do problema, empregando o método dos elementos decontorno (MEC), é apresentado neste trabalho. Devido à presença de atrito e restrições de contato,esse problema torna-se um caso não linear. A não linearidade do problema foi contornada com aaplicação incremental de carga e o uso de um método de resolução de sistemas não lineares. A zonade contato é uma das variáveis do problema e pode conter estados de adesão e escorregamento,simultaneamente. Esses estados dependem dos esforços normais e tangenciais no componente epodem variar durante o processo de aplicação de carga. Dessa forma, cada incremento de cargapode perturbar em relação ao estado anterior. Portanto, o cálculo de variáveis e a atualização dosistema de equações em cada iteração é indispensável. Por este motivo, um algoritmo robusto paradefinição dos estados de contato é proposto. Como o sistema de equações obtido é não linear, o usode um método numérico adequado é exigido. Para a solução deste sistema, o método de Newtonfoi aplicado, o qual permite a verificação do estado de contato em cada incremento. A análiseé feita com o uso de elementos quadráticos contínuos, apresentando resultados contínuos e semoscilação. A comparação dos resultados com as soluções analíticas de Hertz e Mindlin-Cattaneomostram boa concordância.

Palavras-chave: Contato Mecânico, Método dos elementos de contorno, Método de Newton,Elasticidade.

xiii

Abstract

SHATERZADEH-YAZDI, Mohammad Hossein. Contact analysis between two elasticbodies using the Boundary Element Method. 2015. 102p. Thesis (Master’s degree). Faculdade deEngenharia Mecânica, Universidade Estadual de Campinas, Campinas.

The computation of stresses and strains on the components is of great importance, whenthe contact mechanics problems between two elastic bodies are analyzed. In particular cases,bodies are subjected to normal and shear loading in the presence of friction, which increases thecomplexity of the problem. The study of the phenomenon and modeling of the problem, usingthe boundary element method (BEM), are presented in this work. Due to the presence of frictionand natural restrictions, this problem becomes non-linear. The non-linearity of the problem wassolved with an incremental applied load and with the use of solvers to non linear systems. Thecontact zone can contain stick and slip states, simultaneously. These states are dependent on thenormal and shear forces on the component and can vary during the application load process. Thus,each load increment can violate the previous state and therefore, the evaluation of variables andthe updating of the system of equations after each iteration is indispensable. For this reason, arobust algorithm for contact state definition is suggested. Since a non linear system of equationsis obtained, an appropriate numerical method is required. To solve this system, Newton’s methodis applied, which allows the verification of the state of contact at each increment. The analysis isdone with the use of quadratic continuous elements and provides continuous and non-oscillatoryresults. Comparisons of the results with the analytical solutions of Hertz and Mindlin-Cattaneoshow good agreement.

Keywords: Contact Mechanics, The Boundary Element Method, Newton Method, Elasticity.

xv

Lista de Ilustrações

1.1 Modelo de contato fretting com malha MEF (Hojjati Talemi, 2014). . . . . . . . . 62.1 Corpo em equilíbrio sob as forças externas. . . . . . . . . . . . . . . . . . . . . . 102.2 Vetores de tensão sobre as faces do elemento infinitesimal de volume. . . . . . . . 112.3 Placa fina sujeita a um carregamento uniforme distribuído na espessura. . . . . . . 152.4 Cilindro longo sujeito a um carregamento uniforme ao longo do comprimento. . . . 162.5 Teoria da energia de distorção máxima (von Mises). . . . . . . . . . . . . . . . . . 183.1 Caracterização de contatos: (a) Incompleto e não conforme; (b) Completo; (c) In-

completo com singularidade; (d) Incompleto e conforme. . . . . . . . . . . . . . . 203.2 Contato normal de dois corpos elásticos similares . . . . . . . . . . . . . . . . . . 223.3 Semi-plano sujeito a força normal e tangencial. . . . . . . . . . . . . . . . . . . . 243.4 Modo de integração de cargas normal e tangencial no semi-plano. . . . . . . . . . 263.5 Tensões subsuperficiais ao longo do eixo y de simetria. . . . . . . . . . . . . . . . 293.6 Contato de cilindros sujeito uma força normal e uma carga tangencial. . . . . . . . 303.7 Zonas de adesão e escorregamento entre dois cilindros em regime de escorrega-

mento parcial. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.8 Tensões normalizadas de Mindlin-Cattaneo em regime de escorregamento parcial. . 334.1 O ponto fonte x′ pertence ao contorno. . . . . . . . . . . . . . . . . . . . . . . . . 394.2 Ângulo interno no contorno do componente. . . . . . . . . . . . . . . . . . . . . . 414.3 Elementos quadráticos contínuos. . . . . . . . . . . . . . . . . . . . . . . . . . . . 424.4 Coordenada adimensional ξ para elementos quadráticos continuos. . . . . . . . . . 434.5 Funções de forma de elementos quadráticos contínuos. . . . . . . . . . . . . . . . 444.6 Tensões no contorno. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 485.1 Região de contato Γc e região sem contato Γnc. . . . . . . . . . . . . . . . . . . . 515.2 Sistema de coordenadas globais. . . . . . . . . . . . . . . . . . . . . . . . . . . . 605.3 Sistema de coordenadas locais. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 616.1 Ilustração da iteração do método de Newton. . . . . . . . . . . . . . . . . . . . . . 666.2 Demostração de caso de comportamento cíclico. . . . . . . . . . . . . . . . . . . . 677.1 Contato de dois cilindros sob carregamento vertical e horizontal. . . . . . . . . . . 717.2 Modelo de duas sapatas em contato. . . . . . . . . . . . . . . . . . . . . . . . . . 727.3 Modelagem do problema no programa implementado. . . . . . . . . . . . . . . . . 73

xvii

7.4 Tensão normal com 60 elementos quadráticos contínuos na região de contato com-parando com a solução de Hertz no contato sem atrito. . . . . . . . . . . . . . . . . 74

7.5 Tensão normal com com 120 elementos quadráticos contínuos na região de contatocomparando com solução de Hertz no contato sem atrito. . . . . . . . . . . . . . . 74

7.6 Deslocamentos dos corpos ao longo do eixo x e y no contato sem atrito. . . . . . . 757.7 Mapa de cor das tensões equivalente de von Mises nos corpos sob carregamento

vertical sem atrito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 767.8 Variação de tensões principais no plano ao longo do eixo y e comparação com a

solução analítica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 777.9 Tensões normal e cisalhante no corpo 1 com 60 elementos quadráticos contínuos

na região de contato considerando µ = 0, 2. . . . . . . . . . . . . . . . . . . . . . 787.10 Tensões normal e cisalhante no corpo 1 com a malha mais refinada no contato com

atrito . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 797.11 Deslocamentos dos corpos ao longo do eixo x e y no contato com atrito. . . . . . . 807.12 Mapa de cor das tensões equivalentes de von Mises nos corpos sob carregamento

vertical com atrito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 807.13 Modelagem do problema no programa implementado. . . . . . . . . . . . . . . . . 817.14 Condições de contorno na aresta superior. . . . . . . . . . . . . . . . . . . . . . . 827.15 Tensões normal e cisalhante no corpo 1 com a malha menos refinada e comparação

com a solução analítica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 837.16 Tensões normal e cisalhante no corpo 1 com a malha mais refinada e comparação

com a solução analítica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 837.17 Deslocamento dos corpos ao longo do eixo x e y no contato sob carregamento

vertical e horizontal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 847.18 Deslocamento dos corpos na direção y no lado esquerdo do corpo. . . . . . . . . . 857.19 Deslocamento dos corpos na direção y no lado direito do corpo. . . . . . . . . . . 857.20 Mapa de cor das tensões equivalente de von Mises nos corpos sob carregamento

horizontal com atrito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 867.21 Tensões normal e cisalhante no corpo 1 sob carregamento vertical e horizontal e

solução analítica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 877.22 Deslocamento dos corpos ao longo do eixo x e y sob carregamento vertical e hori-

zontal com atrito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 887.23 Mapa de cor das tensões equivalente de von Mises nos corpos sob carregamento

vertical e horizontal com atrito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

xviii

SUMÁRIO

Lista de Ilustrações xvii

SUMÁRIO xix

1 Introdução 11.1 Revisão bibliográfica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2 Motivação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.3 Organização do texto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2 Embasamento Teórico 92.1 Elasticidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.2 Tensão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.3 Deformação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.4 Lei de Hooke . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.5 Estado plano de tensão e deformação . . . . . . . . . . . . . . . . . . . . . . . . . 152.6 Teoria da energia de distorção máxima . . . . . . . . . . . . . . . . . . . . . . . . 16

3 Mecânica do contato 193.1 Introdução ao contato mecânico . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.1.1 Contato sem atrito . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.1.2 Contato com atrito . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.1.3 Modos de contato . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.2 Contato de dois corpos elásticos . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.3 Solução analítica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233.4 Teoria de Hertz para o contato elástico . . . . . . . . . . . . . . . . . . . . . . . . 273.5 Contato de dois cilindros sujeito a escorregamento parcial . . . . . . . . . . . . . . 30

4 O Método dos Elementos de Contorno 354.1 A equação integral de contorno . . . . . . . . . . . . . . . . . . . . . . . . . . . . 354.2 Solução Fundamental . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 374.3 A equação integral de contorno . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

xix

4.4 Elementos de contorno quadráticos contínuos . . . . . . . . . . . . . . . . . . . . 424.4.1 Integração das matrizes H e G quando o ponto fonte não pertence ao elemento 45

4.5 Tensões no contorno . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

5 Método de elementos de contornos para problemas de contato 515.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 515.2 Aplicação incremental de carga . . . . . . . . . . . . . . . . . . . . . . . . . . . . 575.3 Sistema de coordenadas na região de contato . . . . . . . . . . . . . . . . . . . . . 595.4 Restrições na região de contato . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

5.4.1 Adesão ou escorregamento parcial . . . . . . . . . . . . . . . . . . . . . . 625.4.2 Escorregamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 625.4.3 Separação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

5.5 Decisão do estado de contato . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 635.5.1 Separação ou Contato . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 635.5.2 Adesão ou Deslizamento . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

6 Método de solução de sistemas não lineares 656.1 Princípios do método de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . 656.2 Sistema multidimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 686.3 Algoritmo do programa principal . . . . . . . . . . . . . . . . . . . . . . . . . . . 696.4 Algoritmo para método de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . 706.5 Algoritmo para a montagem das matrizes ANL e bNL . . . . . . . . . . . . . . . . 70

7 Simulação numérica e resultados 717.1 Exemplo 1: Contato de duas sapatas elásticas sem atrito sob carregamento vertical . 71

7.1.1 Resultados numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 737.2 Exemplo 2: Contato de duas sapatas elásticas com atrito sob carregamento vertical . 78

7.2.1 Resultados numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 787.3 Exemplo 3: Contato de duas sapatas elásticas sob carregamento vertical e horizontal 81

7.3.1 Resultados numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

8 Considerações Finais 918.1 Conclusões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 918.2 Sugestões para trabalhos futuros . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

Referências 93

xx

APÊNDICES 101

A Código computacional da solução analítica Mindlin-Cattaneo 101

xxi

1 Introdução

A interação mecânica entre corpos em contato tem grande importância nas atividades diárias,especificamente nos problemas de engenharia. Wriggers e Laursen (2006) afirmam: "de modo ge-

ral, podemos dizer todos os movimentos no nosso planeta, envolve contato e atrito". Em uma visãotécnica, o contato entre dois corpos pode ocorrer em um ponto, ao longo de uma linha ou sobre umasuperfície. No caso de carregamentos e configurações complexas, pode-se encontrar uma combi-nação dos contatos mencionados e dependendo da força aplicada, a região de contato pode mudar.Por outro lado, a interação entre duas superfícies de contato é complexa, pois o comportamento decontato é sensível à rigidez do material, à sua textura, ao seu acabamento, à topologia da superfíciede contato, à taxa de carregamento, à magnitude do carregamento, à direção do carregamento emrelação a região de contato, aos suportes do corpo, etc.

Devido à complexidade e importância tecnológica, problemas de contato têm sido muito es-tudados nos últimos anos. Grande parte dos estudos podem ser encontrados nos livros Johnson(1987), Wriggers e Laursen (2006) e Popov (2010). O comportamento na interface de contato émuito influenciado pelo atrito, a transferência de carga e a interação entre duas superfícies de con-tato. O contato é um parâmetro mecânico importante nas estruturas, pois a carga entre os compo-nentes nas máquinas, conjuntos aparafusados, engrenagens, rolamentos, fixação de pás de turbinas,entre outras, é transferida por contato. Calcular o valor de tensões e deformações devido ao contato,especificamente na presença de atrito, é de grande importância prática.

Considerando que a maioria das estruturas metálicas não sobrevive indefinidamente, a análisedo contato mecânico ajuda a evitar falhas catastróficas. Por exemplo, na indústria aeroespacial, ainiciação de trincas nos pontos de concentração de tensão é uma das preocupações dos engenheiros.Os furos dos rebites, por exemplo, são pontos de concentração de tensão sob condições de contato.Em tal situação, a mecânica do contato deve ser utilizada para fornecer as avaliações necessárias eprocedimentos para lidar com esses problemas, tanto na fase de projeto quanto em serviço.

Também podem ser destacadas outras aplicações particularmente úteis: o caso de fretting

que aparece quando o contato mecânico está associado a cargas cíclicas. Esta situação apresentapequenas oscilações dos corpos na região de contato, o que pode acelerar a iniciação de trincassuperficial ou subsuperficial e fazer com que a propagação destas trincas levem os componentes àfalha catastrófica.

1

Em certas situações, uma análise de contato é feita para avaliar a quantidade de dano, comono caso de fretting (Waterhouse, 1981). Em outros casos, o contato é desejado e é estudado paraaperfeiçoar a utilização de componentes mecânicos, tais como as travas mecânicas. A caracteri-zação dos problemas muitas vezes ocorre por meio de equações diferenciais e integrais. No casode contato mecânico com atrito, a natureza não linear do fenômeno torna mais difícil uma solu-ção exata. Neste contexto, é inevitável o uso de um método numérico para resolver o sistema deequações não lineares ou sistemas de equações lineares com restrições de desigualidade.

Atualmente, o método dos elementos finitos (MEF) é amplamente utilizado para realizar aanálise de problemas de contato. O MEF se baseia numa estratégia de aproximação variacionalque discretiza o corpo em elementos de tamanhos finitos. Cada um dos elementos é descrito poruma aproximação por partes das equações governantes utilizando as abordagens variacionais ouresiduais ponderadas. Reunindo as equações para todos os elementos, um sistema simples de equa-ções algébricas que representa o comportamento global do corpo pode ser obtido e resolvido. Odomínio de aplicação bastante amplo do MEF representa um grande desafio para qualquer outrométodo numérico existente. No entanto, uma de suas desvantagens é a necessidade de discretiza-ção de todo o corpo, que pode conduzir a um sistema de matrizes muito grande para ser resolvido,especificamente para problemas tridimensionais complexos.

Apesar de alguns destes problemas terem sido parcialmente resolvidos pelo desenvolvimentorecente de algoritmos rápidos para a solução de sistemas de equações e computadores com grandecapacidade de cálculo, algumas das dificuldades inerentes associadas com o MEF permanecem.Segundo Man (1994), o MEF ainda é ineficiente e demorado para os problemas onde o contornomuda constantemente, tais como aqueles encontrados na mecânica da fratura linear elástica e namecânica do contato devido aos cálculos desnecessários no interior do domínio.

Ao contrário do MEF, o método de elementos de contorno (MEC) evita a discretização detodo o domínio, usando uma abordagem matemática diferente. Esta técnica transforma analitica-mente o conjunto de equações diferenciais lineares governantes em conjunto de equações integraisao longo do contorno do problema. Essa transformação permite usar sistemas de discretizaçãoque envolvem apenas o contorno do corpo como foi dito por Fredholm (1903), Smirnov (1964)e outros. A abordagem do MEC tem muitas vantagens sobre outras técnicas numéricas. Essasvantagens são resumidas como se segue:

2

• Reduz a dimensionalidade do problema em um, resultando em um sistema menor de equa-ções com uma redução considerável nos dados necessários para a análise.

• O MEC considera contínua a modelagem no interior do domínio da solução, uma vez quenenhuma discretização do interior é necessária, o que conduz a uma elevada precisão no cálculodas tensões e deslocamentos interiores.

• As tensões são calculadas com a mesma precisão dos deslocamentos.

• O método é bem adequado para os problemas de domínios infinitos, tais como mecânicados solos, acústica, aeroelasticidade, dentre outros. Nestes casos, os métodos clássicos de domíniosão menos eficientes.

Pode-se argumentar a partir do ponto de vista numérico que, como o contato acontece nocontorno, uma solução de contorno como a do método dos elementos de contorno, em vez de umprocedimento de solução de domínio, seria mais adequada para a análise destes tipos de problemas.Além disso, na análise dos elementos de contorno, tanto os deslocamentos como as tensões sãoobtidas com a mesma precisão. No caso do método dos elementos finitos, por exemplo, as tensõessão obtidas com precisões inferiores aos deslocamentos.

1.1 Revisão bibliográfica

O primeiro modelo para problemas de contato elástico foi proposto por Hertz (1882). Aformulação de Hertz é dada para o caso de contato entre dois corpos elásticos sem atrito, e aindahoje é considerada como uma solução prática para problemas sem atrito (Barbosa e Ghaboussi,1990). Meio século depois, ? resolveu problemas unilaterais (Elástico-rígido).

Cattaneo (1938) publicou a primeira extensão da formulação de Hertz para corpos elásticos,considerando o atrito na superfície de contato de esferas de materiais idênticos sujeitas a desliza-mento tangencial. Uma década depois, Mindlin (1949) estudou o mesmo problema, supondo quetoda região de contato está em adesão. Essa suposição mostrou que, se toda região estiver em estadode adesão, aprecem tensões tangenciais singulares. Para resolver essa singularidade, a aplicação de

3

coeficiente de atrito Coulomb foi sugerida.

Outra extensão da formulação do Hertz foi feita por Goodman (1962). Ele considerou ocontato sem a existência de deslizamento e sugeriu a aplicação gradual de carga normal. Curnier(1984), com uso da analogia de plasticidade, apresentou aspectos da teoria de atrito. No mesmoano, Cocu (1984) demonstrou a junção da teoria de atrito Coulomb com a formulação de Signorini.

Duvaut e Lions (1972) apresentaram uma nova visão sobre o problema que é mais abrangente.Eles formularam o problema como um problema quase estático com atrito na forma de princípiosvariacionais. Como referência para esse princípio, pode-se mencionar o livro do Kikuchi e Oden(1988).

Como exemplo de análise de contato utilizando métodos numéricos, pode-se destacar o tra-balho de Panagiotopoulos (1975) como pioneiro. Entre os trabalhos com o método dos elementosfinitos, pode citar os trabalhos de Singh e Paul (1974), Frangi e Novati (2003) e Batra (1981). Paraanálises feitas com o método dos elementos de contorno, os trabalhos do Andersson (1981), Paris eGarrido (1989), Abascal (1995), Man et al. (1993) e Rodríguez-Tembleque Solano (2009) podemser citados.

No desenvolvimento dos estudos com uso do método lagrangiano aumentado, os trabalhosde Landers e Taylor (1986), Wriggers et al. (1985), Kikuchi e Oden (1988), Alart e Curnier (1991)e Simo e Laursen (1992) são destacados. O contato termomecânico é estudado por Zavarise et al.

(1995) e para problema de fadiga por fretting, Strómberg (1999) foi um dos pioneiros.

Melhorias da formulação lagrangiana com uso de programação matemática foi realizado porKlarbring (1986) e Christensen (1997). Nesses trabalhos, eles destacam vários problemas de con-tato com várias técnicas de programação matemática. Os avanços mais recentes da formulaçãolagrangiana foram feitos com base no método Mortar e MMLs por Puso e Laursen (2004), Rebelet al. (2002) e González et al. (2008).

4

1.2 Motivação

Um dos casos mais específicos, investido e explorado nos últimos anos é caso de fretting.Esse fenômeno se refere ao desgaste e, às vezes, corrosão das superfícies de contato. Estes danossão induzidos sob cargas cíclicas, tais como as induzidas por vibração. O manual da ASM defadiga e fratura, Metals (1961), define fretting como: "Um processo de desgaste que ocorre na área

de contato entre dois materiais sob carga e sujeitos a movimentos relativos devido à vibração ou

alguma outra força."A amplitude do movimento de deslizamento relativo é muitas vezes na ordemde mícron a milímetros, mas pode ser tão baixa quanto 3 a 4 nanômetros.

Muitos materiais de engenharia têm aplicações onde componentes estruturais são submetidosa condições de fretting como, por exemplo, em juntas parafusadas e rebitadas no acoplamento deeixos com engrenagens e/ou rolamentos, na interface da montagem das palhetas com o disco deturbinas ou compressores (Ruiz et al. (1984); Ruiz e Chen (1986); Ruiz e Nowell (2000)), nasjuntas rebitadas da fuselagem de aeronaves (Harish e Farris, 1998), etc. Testes experimentais têmmostrado que a ocorrência da fadiga por fretting pode reduzir em até 90% a resistência à fadiga deum material metálico (McDowell et al., 1954).

Uma das abordagens de estudo de fadiga por fretting considera que o fenômeno poderiaser tratado como um problema de fadiga convencional na presença de um concentrador de ten-são (notch analogue). Com isto, minimiza-se a consideração do efeito do desgaste superficiale maximiza-se o efeito de concentração de tensões na região do contato. Giannakopoulos et al.

(2000) observaram que o campo de tensão resultante do contato entre uma sapata plana com cantosarredondados e um semi-plano (Figura 1.1) é similar ao campo de tensão de corpos entalhados esugeriram que se deveria explorar esta característica para estabelecer metodologias de previsão devida ou resistência à fadiga por fretting.

Deve ser lembrado que para modelar esse tipo de problema, requer-se consideráveis cuida-dos: não somente refinar a malha cuidadosamente nas vizinhanças do contato, mas também muitosoutros aspectos devem ser analisados como convergência dos resultados. O não refinamento ade-quado da malha pode inserir certas imprecisões na solução, particularmente na posição das regiõesde escorregamento e adesão do contato, as quais são relativamente difíceis de localizar. Sendo as-sim, para estudos com fins de caracterização fenomenológica da fadiga sob condições de fretting, épreferível utilizar dados de testes que empreguem geometrias idealizadas e bem definidas, de modo

5

Figura 1.1: Modelo de contato fretting com malha MEF (Hojjati Talemi, 2014).

que a natureza do contato e das tensões/deslocamentos induzidos pelo contato seja facilmente con-trolável, possuam repetibilidade e pequena sensibilidade às imperfeições de fabricação.

Muitos fatores que influenciam a resistência à fadiga por fretting, como a pressão no contato,a amplitude do escorregamento relativo, condições ambientais e materiais, ainda não tinham sidoavaliados completamente. Entretanto, em 1968, surgiram os trabalhos de Nishioka et al. (1968),seguidos pelas publicações de Nishioka e Hirakawa (1969) e Nishioka e Kenji (1972), que exami-naram a influência desses fatores de maneira independente. Uma das principais conclusões destesestudos foi que havia uma faixa de deslocamentos tangenciais que acelerava o processo de fadigapor fretting.

Bramhall (1973) observou o efeito do tamanho do contato na vida à fadiga, após a realizaçãode uma série de experimentos onde se mantinha o estado de tensão superficial constante de testepara teste, mas variava-se o tamanho do contato. Para qualquer tamanho de contato inferior a umtamanho crítico observou-se que a vida era infinita (> 107ciclos), enquanto que para maiores ta-manhos de contato, a falha ocorria. Posteriormente, outros pesquisadores como Nowell (1988) eAraujo (2000), confirmaram a existência deste efeito para outros materiais.

Fouvry et al. (1998) utilizaram de experimentos com contatos esfera-plano sob condiçõesde escorregamento parcial para validar a aplicação de alguns critérios de fadiga multiaxial e veri-ficaram que os resultados obtidos não eram satisfatórios quando o campo de tensões apresentavaseveros gradientes de tensão. Araujo e Nowell (2002) conduziram uma abordagem similar e ve-

6

rificaram que melhores resultados poderiam ser obtidos utilizando uma zona de processo que nãopareceu a princípio ser característica própria do material.

Soluções numéricas para problemas de contato estão disponíveis, geralmente usando o mé-todo dos elementos finitos (MEF) como, por exemplo, nos trabalhos de Johansson (1994), Stróm-berg et al. (1996), Omberg (1997) e Strómberg (1999). O MEF também foi utilizado por Iremanet al. (2009) que apresentaram as formulações e os fundamentos para simulações de desgaste de-vido ao fretting. Madge et al. (2007) propuseram um modelo 2D dos elementos finitos para afadiga por fretting. Mais recentemente, pode-se mencionar a análise numérica e experimental deproblemas de desgaste por fretting apresentada por Wei et al. (2011), Páczelt et al. (2012), Bai-etto et al. (2013), Giner et al. (2014), (Hojjati Talemi, 2014). Também podem ser citados algunspoucos trabalhos no qual se usa o Método dos Elementos de Contorno (MEC) como, por exemplo,Rodríguez-Tembleque et al. (2012).

Assim, neste trabalho será estudado o fenômeno de contato de dois cilindros como primeiraetapa de modelagem de contato fretting. Este modelo é utilizado para calibrar os modelos numéri-cos. Neste caso, o campo de tensão possui solução analítica bem definida (Hills e Nowell, 1994) efoi adotada por outros pesquisadores (Nowell (1988); Araujo (2000); Fouvry et al. (2002)).

1.3 Organização do texto

A organização deste trabalho se encontra da seguinte forma:

Capítulo 1: Apresenta uma introdução sobre problemas de contato, motivação do estudo doproblema, objetivos e uma revisão bibliográfica sobre o caso.

Capítulo 2: Tratando de um problema de elasticidade, esse capítulo apresenta os concei-tos básicos da área de mecânica de contínuo que foram utilizados no desenvolvimento dotrabalho.

Capítulo 3: A formulação do método dos elementos de contorno com foco nos elementosquadráticos é apresentada neste capítulo.

Capítulo 4: A aplicação do MEC para problemas de contato, construção de matrizes e sis-temas, restrições e pontos importantes a serem considerados na implementação do problema

7

se encontram nesse capítulo.

Capítulo 5: Este capítulo apresenta uma breve introdução sobre a solução de sistemas nãolineares e também a formulação do método empregado.

Capítulo 6: Os exemplos modelados e respectivos resultados são apresentados e analisados.

Capítulo 7: Conclusões, últimos comentários e trabalhos futuros são apresentados.

Anexos: O código da solução analítica de Mindlin-Cattaneo é apresentado.

8

2 Embasamento Teórico

2.1 Elasticidade

Todos os materiais estruturais possuem, em certa medida, um módulo de elasticidade. Desdeque as forças externas, que produzem uma deformação na estrutura, não excedam um certo limite,as deformações desaparecem com a remoção das forças, pois as tensões são proporcionais às de-formações na fase elástica do material (Timoshenko e Goodier, 1951). Neste estudo será assumidoque os corpos submetidos à ação de forças externas são perfeitamente elásticos e retomam a suaforma inicial completamente após a remoção das forças.

A estrutura molecular dos corpos elásticos não será considerada aqui. Será assumido que amatéria de um corpo elástico é homogênea e distribuída continuamente ao longo do seu volume,de modo que um pequeno elemento do corpo possui as mesmas propriedades físicas do corpo. Parasimplificar a discussão, também será assumido que o corpo é isotrópico, ou seja, as propriedadeselásticas são as mesmas em todas as direções.

2.2 Tensão

Tensão é uma medida das forças internas que atuam em pontos dentro de um corpo. Quantita-tivamente, é uma medida da força média por unidade de área de uma superfície no interior do corposobre a qual atuam as forças internas. Estas forças internas surgem como reação às forças externasaplicadas ao corpo. Uma vez que o corpo deformável carregado é assumido como contínuo, estasforças internas são distribuídas de forma contínua dentro do volume do corpo e o corpo tem umadeformação contínua.

A Figura 2.1 representa um corpo em equilíbrio sob as forças P1,...,P7. Imagine o corpodividido em duas partes A e B através de um corte na seção mm. Será assumido que estas forçasestão distribuídas sobre a área mm continuamente, da mesma forma que a pressão é continuamentedistribuída sobre a superfície sob a qual atua.

A área ∆S é uma área infinitesimal com normal unitária ni e ∆F é uma força interna exercida

9

A

B

O

N

z

x

y

m

mP1

P2

P3

P4

P5

P6

P7

Figura 2.1: Corpo em equilíbrio sob as forças externas.

sobre essa área. Define-se o vetor tensão como:

Ti(n) = lim

∆S→ 0

∆Fi∆S

=dFidS

(2.1)

O vetor de tensão Ti depende da sua localização no corpo e da orientação do plano sobre oqual atua. Dependendo da orientação do plano considerado, o vetor de tensão pode não ser neces-sariamente perpendicular a esse plano e pode ser decomposto em dois componentes:

• Normal ao plano, chamado de tensão normal:

σn = lim∆S→ 0

∆Fn∆S

(2.2)

onde ∆Fn é o componente normal da força ∆F que atua sobre a área ∆S.

• Paralelo a este plano, chamado tensão de cisalhamento:

τ = lim∆S→ 0

∆Fs∆S

(2.3)

onde ∆Fs é a componente tangencial da força ∆F que atua sobre a área ∆S. A tensão de cisalha-mento pode ser ainda decomposta em dois vetores mutuamente perpendiculares.

10

Em um elemento infinitesimal de volume do corpo elástico, o vetor de tensão pode ser escritocomo três vetores de tensão em cada face nas direções de coordenadas cartesianas (x, y, z). Umvetor normal à superfície e outros dois tangentes à superfície e perpendiculares entre si. Portanto,observam-se 9 vetores de tensão em um volume de corpo.

σz

τzxτzy

σy

τyx

τyz

σx τxy

τxz

x

z

y

Figura 2.2: Vetores de tensão sobre as faces do elemento infinitesimal de volume.

As tensões sobre um volume infinitesimal de um corpo podem ser organizadas em uma matrizsimétrica pois há equilíbrio de momentos. Utilizando-se notação indicial, pode-se escrever cadatensão como σij , onde i representa a coordenada normal ao plano em que a tensão se encontra ej é a direção do componente do vetor tensão. Essa matriz é denominada de tensor de tensões deCauchy e é representada pelo símbolo σ, como abaixo:

σ =

σxx σxy σxz

σyx σyy σyz

σzx σzy σzz

=

σx τxy τxz

τyx σy τyz

τzx τzy σz

(2.4)

O teorema de Cauchy afirma que existe um campo tensorial σ, chamado de tensor de Cauchy,independente de n, tal que T é uma função linear de n:

T (n) = σn (2.5)

11

Definindo bi como uma força por unidade de volume, tem-se:

bi =

bxbybz

(2.6)

Portanto o equilíbrio estático pode ser escrito como:

σij,j + bi = 0 (2.7)

2.3 Deformação

A deformação de um corpo elástico é definida como a variação de comprimento e forma emcerta direção. Se considerar u como a variação de comprimento e dividir pelo comprimento inicialobtém-se a deformação linear média. Em um elemento infinitesimal, a deformação é definida comoa variação de comprimento u por unidade de comprimento, e pode ser escrita como:

εxx =∂ux∂x

, εyy =∂uy∂y

, εzz =∂uz∂z

(2.8)

As Equações (2.8) representam as deformações normais ou lineares. As componentes ci-salhantes de deformação também fazem parte da deformação. Estas componentes são chamadasdeformações angulares e podem ser escritas como:

εxy =1

2(∂ux∂y

+∂uy∂x

), εyz =1

2(∂uy∂z

+∂uz∂y

) e εzx =1

2(∂uz∂x

+∂ux∂z

) (2.9)

Todos os componentes de deformação podem ser escritos na forma de uma matriz denomi-nada tensor de deformações representado por ε que envolve todos os componentes normais e decisalhamento do tensor de deformação.

12

ε =

εxx εxy εxz

εyx εyy εyz

εzx εzy εzz

(2.10)

Na matriz (2.10), devido ao equilíbrio de momentos, os componentes de cisalhamento dadeformação são iguais em dupla, por simetria, ou seja εij = εji.

As equações (2.8) e (2.9) podem ser reescritas como:

ε(u) = ∇su =1

2(∇u+∇uT ) (2.11)

2.4 Lei de Hooke

Como já foi explicado anteriormente, qualquer material, sobre o qual atua uma força, sofreráuma deformação. Em materiais com comportamento elástico linear, a deformação do corpo serelaciona diretamente com as tensões atuantes no corpo (Hibbeler, 2004). Esta relação é dada pelalei de Hooke que relaciona o tensor de tensões de Cauchy e o tensor de deformações.

No caso de material homogêneo e isotrópico, a lei de Hooke pode ser escrita de forma gene-ralizada como:

σij = λδijuk,k +G(ui,j + uj,i) (2.12)

onde λ e G são constantes de Lamé, expressas em termos do módulo Young E e do coeficiente dePoisson ν e são definidos assim:

λ =νE

(1 + ν)(1− 2ν)(2.13)

13

G =E

2(1 + ν)(2.14)

G é chamado de módulo de cisalhamento e δij é o delta de Kronecker cujas propriedades são:

δij =

0 i 6= j

1 i = j(2.15)

e uk,k é o traço do tensor de deformação em notação indicial, ou seja:

uk,k = εkk = ux,x + uy,y + uz,z (2.16)

A inversa da Equação (2.12) pode ser escrita como:

εij =1

2G[σij −

ν

(1 + ν)σkkδij] (2.17)

As equações (2.12) e (2.17) também podem ser reescritas em termos de E e ν como a seguir:

σij =E

(1 + ν)[εij +

ν

(1− 2ν)εkkδij] (2.18)

eεij =

(1 + ν)

E[σij +

ν

(1 + ν)σkkδij] (2.19)

Derivando a Eq. (2.12), substituindo na equação de equilíbrio (2.7) e substituindo os compo-nentes de deformação com as derivadas de deslocamentos obtidos a partir da Eq.(2.9), obtém-se aequação de Navier para equilíbrio estático, dada por:

uj,ii +1

(1− 2ν)ui,ij +

1

Gbi = 0 (2.20)

Esta equação é particularmente conveniente se as condições de contorno de deslocamento sãoespecificadas. Da mesma forma, as condições de contorno de forças de superfície podem ser dadaspelo vetor de tensões definido como:

14

Ti = σijnj = (λδijuk,k)nj + 2Gεij.nj (2.21)

onde nj indica o cosseno diretor do vetor normal apontando para o exterior do domínio do corpo.

2.5 Estado plano de tensão e deformação

Estado plano de tensão: Este modo é considerado quando a espessura do componente érelativamente fina comparando com outras dimensões, e a força é uniformemente distribuídaao longo da espessura, como mostrado na Figura 2.3.

yy

zx

Figura 2.3: Placa fina sujeita a um carregamento uniforme distribuído na espessura.

No estado plano de tensão as componentes σz, τzz e τzy são zero em ambos os lados docorpo e também ao longo das faces. Portanto, é suficiente calcular as tensões σx, σy e τxy docomponente.

Estado plano de deformação: Outra simplificação possível é quando a dimensão do compo-nente na direção z é muito grande e tem seção uniforme. Nesse caso, o carregamento aplicadoestá distribuído de forma uniforme ao longo do eixo z do corpo. Assim sendo, será assumidoque a seção inteira tem a mesma condição. Portanto, nos casos como por exemplo de umcilindro longo (Figura 2.4), considera-se deformação ao longo do eixo z igual zero.

Desta forma, basta analisar uma fatia do corpo ao longo do eixo z. Os componentes u e v dedeslocamento são funções de x e y, sendo independentes da coordenada z. Uma vez que o

15

y

x

y

z

Figura 2.4: Cilindro longo sujeito a um carregamento uniforme ao longo do comprimento.

deslocamento w é zero, obtém-se:

γyz =∂v

∂z+∂w

∂z= 0

γzz =∂u

∂z+∂w

∂z= 0

εx =∂w

∂z= 0 (2.22)

A Equação (2.22) define que no estado plano de deformação as componentes τxz e τyz datensão são zero. Considerando εz = 0 por meio da lei de Hooke, a tensão σz pode ser escritaem termos de σx e σy da seguinte forma:

σz − ν(σx + σy) = 0 ⇒ σz = ν(σx + σy) (2.23)

Logo, conclui-se que no estado plano de deformação, basta calcular as tensões σx, σy e τxy,que são funções de x e y e dessa forma a tensão σz vai ser uma função de σx e σy.

2.6 Teoria da energia de distorção máxima

Os componentes mecânicos e os elementos estruturais são projetados e fabricados para resis-tir às condições de carregamento que é definido como tensão de escoamento para materiaisdúctil. A tensão de escoamento pode ser obtida num ensaio de tração do corpo de prova sobcarregamento uniaxial. Contudo, quando o componente estiver sujeito a tensão multiaxial, aavaliação da resistência do componente será mais difícil.

16

Existem várias teorias para prever a falha de um componente sujeito a tensão multaxial, comoteoria de tensão de cisalhamento máxima e a teoria da densidade de energia de distorção má-xima. Para usar essas teorias, as tensões principais no ponto desejado devem ser encontradas.Nessa seção, somente a teoria da energia de distorção, também conhecida como critério devon Mises, será abordada pois somente esta teoria foi utilizada no trabalho.

Matemático von Mises (1913) propôs que o escoamento dos matérias começa quando o se-gundo invariante de tensão J2 atinge um valor crítico conhecido como tensão de escoamento.O segundo invariante de tensão, J2, também conhecido como J2 plasticidade ou teoria deescoamento J2 é definido da seguinte forma:

J2 =1

2sijsji = −(s1s2 + s2s3 + s3s1) =

1

6

[(σ1 − σ2)2 + (σ2 − σ3)2 + (σ3 − σ1)2] (2.24)

Matematicamente a função de von Mises é definido como:

f(J2) = J2 − k2 = 0

onde k é limite de escoamento do material sob tensão de cisalhamento puro. O valor de k é√3 vezes menor que limite de escoamento do material, ou seja:

k =σe√

3

Substituindo J2 no tensor tensão de Cauchy, tem-se:

σ2e =

1

2

[(σ1 − σ2)2 + (σ2 − σ3)2 + (σ3 − σ1)2 + 6

(σ2

23 + σ231 + σ2

12

)](2.25)

Quando um corpo está sujeito a uma carga externa, ele se deforma e armazena energia emtodo o volume. Essa energia interna, por unidade de seu volume, é denominada densidade deenergia de deformação e é definida da seguinte forma:

u =1

2σ1ε1 +

1

2σ2ε2 +

1

2σ3ε4 (2.26)

Simplificando a Equação (2.26) para caso o de material com comportamento linear elásticoe utilizando a lei de Hooke, para o estado plano de tensões obtém-se:

u =1

2E

[σ2

1 + σ22 + σ2

3 − 2υ (σ1σ2 + σ1σ3 + σ3σ2)]

(2.27)

17

A parte da densidade de energia que pode provocar a distorção (mudança de forma) doelemento é a porção remanescente das tensões principais com a tensão principal média(σi−σmd) chamada de densidade de energia de distorção. Fazendo esta subtração na Equação(2.27) e simplificando, obtém-se:

ud =1 + υ

3E

(σ2

1 − σ1σ2 + σ22

)(2.28)

Considerando σ1 = σe e σ2 = σ3 = 0, a densidade de energia de distorção para um estadouniaxial de tensão é dada por:

ud =1 + υ

3Eσ2e (2.29)

Pela teoria tem ud = ud no momento do escoamento e, portanto, no caso de tensão plana, ocritério de von Mises tem a forma seguinte:

(σ2

1 − σ1σ2 + σ22

)= σ2

e (2.30)

A Equação (2.30) do critério von Mises é a equação de uma curva elíptica, como é mostradana Figura 2.5. Essa curva define a fronteira entre pontos seguros (pontos dentro da elipse) epontos que causam escoamento no elemento (pontos fora da elípse).

σ1

σ2

σe

σe

-σe

-σe

Figura 2.5: Teoria da energia de distorção máxima (von Mises).

18

3 Mecânica do contato

3.1 Introdução ao contato mecânico

Mecânica do contato é o estudo de deformações de corpos sólidos que se tocam em umponto, ao longo de uma linha ou ao longo de uma superfície. Os aspectos físicos e matemáticosdo contato são desenvolvidos com base na mecânica dos materiais e mecânica do contínuo paramateriais com comportamento elástico, viscoelástico e plástico em contato estático ou dinâmico. Omotivo principal nos estudos de contato sob carga normal é calcular a pressão e a adesão atuandonos corpos na direção perpendicular à superfície e as forças tangenciais de superfície gerada peloatrito entre os componentes.

A primeira etapa para estudar problemas de contato é caracterizar o contato e conhecer emqual categoria ela se enquadra. A classificação a seguir automaticamente oferece informações so-bre a natureza de contato. O primeiro caso (Figura 3.1 (a)) é o contato de um cilindro levementepressionado em um semi-plano elástico. No primeiro instante, os corpos se tocam ao longo de umalinha e conforme se aumenta a carga aplicada, a superfície de contato aumenta. Esse tipo de con-tato se chama contato incompleto, onde a superfície de contato não é conhecida a priori e dependeda força aplicada. A aproximação para o semi-espaço é valida para calcular tensões e deformaçõesuma vez que a meia-largura do contato, conhecido como a, é muito menor que o raio,R, do cilindro(a << R). No segundo caso, a Figura 3.1(b) mostra o contato entre um bloco retangular com umsemi-espaço elástico que é conhecido como contato completo, pois o tamanho da região de contatonão depende da carga aplicada. No contato completo, as bordas do componente não são contínuase por isso a pressão nas bordas é singular.

No terceiro caso (Figura 3.1(c)), o componente tem uma borda abrupta e uma borda comcurvatura suave. Nesse caso, a pressão é singular na borda abrupta e na outra borda tende para zerogradualmente. Além disso, o tamanho da região de contato depende da carga aplicada. O últimocaso (Figura 3.1(d)) é um caso de contato incompleto, onde um cilindro é pressionado contra umfuro com diâmetro um pouco maior que diâmetro do cilindro. Como a largura do contato é umafração considerável do raio do cilindro, o furo não pode ser mais aproximado por um semi-plano.Para resolver esse tipo de problema, deve-se usar formulação apropriada para disco e um plano infi-nito contendo o furo que se torna um caso mais complexo. Essas formulações são apresentadas nos

19

P P P P

(a) (b) (c) (d)

Figura 3.1: Caracterização de contatos: (a) Incompleto e não conforme; (b) Completo; (c) Incom-pleto com singularidade; (d) Incompleto e conforme.

trabalhos do Persson (1964) e Mostofi e Gohar (1980). Pela investigação experimental de Fessler eOllerton (1957), a aproximação por semi-plano é valida para taxa de a/R até 0.3.

Além da classificação apresentada, pelo aspecto das forças tangencias entre os corpos, osproblemas de contato se classificam em dois tipos: problemas de contato sem atrito e problemas decontato com atrito, explicados detalhadamente a seguir.

3.1.1 Contato sem atrito

O contato sem atrito é um contato idealizado que tem aplicação muito restrita. De modo geral,peças bem lubrificadas podem ser modeladas sem atrito.

Na situação de contato sem atrito, os corpos em contato podem deslizar sem resistência aolongo da direção tangencial (paralela à superfície de contato). Devido à carga aplicada e a ausênciado atrito, existe somente força normal de compressão na região de contato. Neste caso, os corpospodem se separar, mas não vão interpenetrar. Nessa situação, as tensões na direção tangencialsempre são nulas e a continuidade da tensão normal dentro da zona de contato é sempre preservada.Além disso, problemas de contato sem atrito são independentes da história do carregamento.

20

3.1.2 Contato com atrito

O atrito é um fenômeno físico encontrado naturalmente em problemas de contato reais.Quando o atrito é levado em consideração, o problema se torna mais complicado. O atrito in-fluencia significativamente o comportamento na região de contato. Por exemplo, o movimento dedeslizamento na direção tangencial de um ponto de contato será limitado pelas forças de atrito tan-genciais (cisalhamento), no ponto de contato, que por sua vez dependem da componente normaldas forças (tensões normais) exercidas no mesmo ponto. A relação entre as componentes tangenci-ais e normais das forças impõe um comportamento não linear entre o movimento de deslizamentodas superfícies em contato e a carga externa.

Em situação de atrito, as condições de contato ou são de adesão (sem deslocamento relativona direção tangencial) ou de escorregamento (com deslizamento na direção tangencial).

3.1.3 Modos de contato

A região em que os contornos podem entrar em contato é chamada de "área potencial decontato". O tamanho desta área depende do problema envolvido, uma vez que é determinadapela geometria do problema e pela magnitude da carga final aplicada. A situação de separaçãoou contato é descrita pelos modos de contato. Segundo Man (1994), os modos de contato em umponto se classificam em:

1- Modo de separação: É definido quando os corpos permanecem separados.

2- Modo de escorregamento: É definido quando os pares de nós não estão restritos na dire-ção tangencial, mas livre para deslizar um sobre o outro.

3- Modo de adesão: É definido quando os pares de nós estão restritos na direção normal etangencial, ou seja, não tem qualquer deslizamento durante uma dada etapa do carregamento.

4- Modo Misto: Os modos de adesão e escorregamento podem ocorrer simultaneamente emuma dada região. Neste caso, o modo é chamado de modo misto.

21

3.2 Contato de dois corpos elásticos

Considere dois corpos elásticos similares com o mesmo módulo de elasticidade sob contatonormal (Figura 3.2 ). O contato entre os corpos desenvolve a pressão de contato em ambos oscorpos e as superfícies são comprimidas e sofrem deslocamento paralelo à superfície livre.

Corpo 1

Corpo 2

P

P

y

x

Figura 3.2: Contato normal de dois corpos elásticos similares

No entanto, como os corpos são elasticamente similares, os pontos relativos dos dois corpossofrerão mesmo deslocamento e não haverá tendência de escorregamento relativo. Se tiver forçastangencias suficientes para causar escorregamento, forças de superfícies limitadas ao coeficiente deatrito surgirão. A relação entre forças tangenciais e forças normais é apresentada na Equação (3.1):

|q(x, y)| = −µp(x, y) (3.1)

onde q é tensão cisalhante, p é pressão normal (p < 0) e µ é a coeficiente de atrito. A tensãocisalhante causa deslocamentos tangencias em cada superfície. Uma vez que os corpos têm propri-edades similares, os deslocamentos normais dos corpos são iguais e a distribuição de pressão nãose altera. Em alguns casos de contato, como por exemplo fretting, a força tangencial não é grandesuficiente para causar escorregamento total de um corpo sobre o outro. Nesses casos, a região decontato é formada de uma região de adesão, onde os dois corpos são acoplados um no outro, e re-gião de escorregamento, onde ocorrem movimentos tangenciais relativos e existe tensão cisalhante

22

relacionado com coeficiente de atrito.

Na zona de escorregamento, a tensão cisalhante é dada pela Equação (3.1) e a direção delaé oposta à direção de deslocamento tangencial no último incremento. Em razão disso, a tensãocisalhante na zona de adesão sempre é menor que o valor relacionado com µ, ou seja:

|q(x, y)| < −µp(x, y) (3.2)

Vale ressaltar novamente que as tensões cisalhantes dos corpos tem a mesma magnitude mascom direções opostas. Nos problemas de contato, a aplicação de carga deve ser de forma iterativa,uma vez que as zonas de adesão e escorregamento são funções da carga aplicada. As zonas devemser recalculadas em cada iteração. Em outras palavras, tanto as zonas de adesão e escorregamentocomo a região efetiva de contato podem mudar de uma iteração para outra.

3.3 Solução analítica

O primeiro passo para resolver problemas de contato é resolver equações integrais e calcular adistribuição de tensões. Suponha o semi-plano da Figura 3.3 sujeito de carga normal P e tangencialQ por unidade de comprimento.

Considerando o semi-plano no estado plano de deformação, as tensões são dadas pela equa-ção de tensão de Airy na forma seguinte:

φ(r, θ) = −rθπ

(Psenθ +Qcosθ) (3.3)

Substituindo a Equação (3.3) na equação biharmonica (Timoshenko e Goodier, 1951), astensões são obtidas da maneira seguinte:

23

P

Q

r

θ

σ

y

x

Figura 3.3: Semi-plano sujeito a força normal e tangencial.

σrr = − 2

πr(Pcosθ −Qsenθ)

σθθ = τrθ = 0 (3.4)

Substituindo a Equação (3.4) na lei de Hooke, as deformações no plano são obtidas como:

εrr =1

8µσrr (1 + κ) + σθθ (κ− 3)

εθθ =1

8µσθθ (1 + κ) + σrr (κ− 3)

γrθ =τrθ/G (3.5)

onde κ = 3 − 4ν para estado plano de deformação, G é módulo de cisalhamento e ν é coeficientede Poisson. Uma vez obtida as deformações do componente, os deslocamentos podem ser obtidose tem a seguinte forma:

24

2µur =− P

2π(κ− 1)θsenθ − cosθ + (κ+ 1) ln rcosθ (3.6)

− Q

2π(κ− 1)θcosθ + senθ − (κ+ 1) ln rsenθ+ C1 (3.7)

2µuθ =P

2π(κ− 1)θcosθ − senθ − (κ+ 1) ln rsenθ (3.8)

+Q

2π−(κ− 1)θsenθ + cosθ − (κ+ 1) ln rcosθ+ C2 (3.9)

onde C1 e C2 são constantes arbitrários cujo valor depende da especificação do deslocamento .Considerando θ = ±π/2, e convertendo para coordenadas cartesianas, para os deslocamentos nasuperfície tem-se:

u = −P(κ− 1

8µ

)sgn(x) +Q

(κ+ 1

4πµ

)ln |x|+ C1

2µ(3.10)

v = −P(κ+ 1

4πµ

)ln |x|+Q

(κ− 1

8µ

)sgn(x) +

C2

2µ(3.11)

As Equações (3.10) e (3.11) são as bases para a solução de problemas de contato e obtençãodos deslocamentos normais e tangencias dos corpos. Todavia, os constantes C1 e C2 são inconve-nientes e geralmente é preferível uso de derivadas de deslocamentos. Derivando as equações 3.10e 3.11 e simplificando, obtém-se:

1

A

∂g

∂x=

1

π

∫ a

−a

q(ξ)dξ

ξ − x+ βp(x) (3.12)

1

A

∂h

∂x=

1

π

∫ a

−a

p(ξ)dξ

ξ − x|x| ≤ a (3.13)

onde h = v1(x)− v2(x), g = u1(x)−u2(x), ξ é a distância do ponto de referência (Figura 3.4), a émetade do comprimento da região efetiva de contato e A para estado plano de deformação é obtida

25

pela Equação (3.14):

A = 2

1− ν2

1

E1

+1− ν2

2

E2

(3.14)

ry

x σ

ξ

dξ

p(ξ)

q(ξ)

a a

Figura 3.4: Modo de integração de cargas normal e tangencial no semi-plano.

Para obter a solução de p(x) na função de h(x), basta inverter a Equação (3.13) e obter aEquação (3.15). Essa equação foi demonstrada por Muskhelishvili (1977):

p(x) = −w(x)

Aπ

∫ a

−a

h′(ξ)dξ

w(ξ)(ξ − x)− Cw(x) (3.15)

onde h′(x) = ∂h/∂x e a solução fundamental w(x) depende da natureza de contato. Para contatode cilindro (não singular) o valor de w(x) é dado por:

w(x) =√a2 − x2 (3.16)

26

3.4 Teoria de Hertz para o contato elástico

A distribuição de pressão, deformação e a área de contato entre dois corpos elásticos circula-res foi estudada por Hertz e é conhecida como teoria de Hertz (1882). Essa teoria é válida para ocontato de corpos não conformes onde a região de contato é muito menor que o tamanho e o raiode curvatura dos corpos (Johnson, 1987). Um caso de aplicação de Hertz é no contato entre doiscilindros longos com eixos paralelos e sem atrito nas superfícies. Nesse caso, o problema pode seraproximado por um modelo plano e resolvido com a teoria de Hertz. Considerando o raio de cur-vatura do corpo como R e a região de contato como a, as suposições necessárias para a aplicaçãoda teoria de Hertz são as seguintes:

As superfícies são contínuas e não conformes (a << R);

A deformação é pequena (não é aplicável para materiais com módulo de elasticidade baixo);

Cada componente pode ser modelado como um semi-espaço elástico (a << R);

As superfícies são sem atrito (qx = qy = 0).

Segundo Hills e Nowell (1994), uma vez respeitadas as condições acima, a aproximação doscorpos pode ser calculada através da Equação (3.17):

h(x) = ∆− 1

2kx2 (3.17)

onde k é a curvatura relativa dos corpos e é dada por:

k =1

R1

+1

R2

(3.18)

ondeR1 eR2 são os raios dos corpos em contato. Como foi dito na seção anterior, apenas a derivadade deslocamento é utilizada, ou seja:

dh

dx= −kx (3.19)

Substituindo as Equações (3.19) e (3.16) na Equação (3.15), obtém-se:

27

p(x) = −√a2 − x2

Aπ

∫ a

−a

kξdξ√a2 − ξ2(ξ − x)

(3.20)

Resolvendo a integral (3.20), conclui-se que, devido à força normal estática, a distribuição depressão é de forma elíptica:

p(x) = − kA

√a2 − x2 (3.21)

Contudo, a meia largura de contato a ainda é desconhecida. Ela pode ser obtida fazendo oequilíbrio entre a pressão de contato e a força aplicada P , ou seja:

P = −∫ a

−ap(ξ)dξ =

πka2

2A(3.22)

Portanto, a meia largura de contato é obtida pela seguinte fórmula:

a2 =2PA

πk(3.23)

E a Equação (3.21) pode ser reescrita como:

p(a) = −p0

√1− (x/a)2 (3.24)

onde p0 é a pressão de contato máxima dada pela equação:

p0 =2P

πa(3.25)

Uma vez obtida à pressão distribuída no corpo (Equação 3.24), as tensões podem ser calcu-ladas. Na interface de contato a tensão total é σx = −p(x) e fora da região de contato é nula. Avariação de tensões principais ao longo do eixo y é dada por:

28

σx = −p0

a

(a2 + 2y2

)√(a2 + y2)− 2y

(3.26)

σy = −p0

a

√(a2 + y2) (3.27)

E consequentemente para tensão cisalhante principal, tem-se:

τ = −p0

a

y − y2

√(a2 + y2)

(3.28)

A Figura 3.5 mostra as curvas de tensões principais e a tensão cisalhante principalao longodo eixo de simetria no meio do componente.

Figura 3.5: Tensões subsuperficiais ao longo do eixo y de simetria.

29

3.5 Contato de dois cilindros sujeito a escorregamento parcial

Na seção anterior o contato sem atrito sob carga vertical foi discutido e as formulações foramdadas. Nesta seção, será analisado a aplicação subsequente ds carga vertical Q menor que a limita-ção por atrito (Q < µP ), conforme na Figura 3.6. Vale reforçar que considerações necessárias parateoria de Hertz continuam necessárias para caso de escorregamento parcial.

P Q

QP

Figura 3.6: Contato de cilindros sujeito uma força normal e uma carga tangencial.

Como foi explicado anteriormente, quando se tem a presença de carga tangencial, a região decontato pode ser dividida em zona de adesão e zona de escorregamento. A divisão entre as zonasnão é conhecida inicialmente e portanto, é necessário definir cada zona.

Inicialmente assume-se a ausência de escorregamento antes da aplicação da força tangencial.Consequentemente, a tensão cisalhante é zero e não haverá movimento relativo entre os corpos, ouseja, toda a região de contato x ≤ a será em modo de adesão e a tensão cisalhante é definida como:

1

π

∫ a

−a

q(ξ)

ξ − xdξ = 0 |x| ≤ a (3.29)

Resolvendo a Equação (3.29), obtém-se:

q(x) =Q

π√a2 − x2

(3.30)

30

Pode-se observar quando x→ ±a, a distribuição de pressão, p(x), dada pela Equação (3.24)tende para zero e a tensão cisalhante, q(x), tende para infinito. Por isso, nas extremidades da zonade contato, a taxa de q(x)/p(x) vai para infinito e somente um coeficiente de atrito infinito podeimpedir o escorregamento do corpo. Em outras palavras, os pontos do corpo nas extremidades daregião de contato são sujeitos a uma tensão singular.

A existência de tensão cisalhante singular em ambas as extremidades sugere a ocorrência deescorregamento. No caso de contato de cilindros onde os corpos têm simetria de geometria e decarga aplicada, é sugerida uma região de adesão no meio da região de contato cercado por duasregiões de escorregamento simétricas. A região de adesão é definida como x ≤ c e a região deescorregamento é onde c < |x| < a, conforme ilustrado na Figura 3.7.

P Q

-c c a-a

Escorregamento Adesão Escorregamento

Figura 3.7: Zonas de adesão e escorregamento entre dois cilindros em regime de escorregamentoparcial.

Para calcular a tensão cisalhante, acrescenta-se uma variável de q′(x) na Equação (3.30) daseguinte forma:

q(x) = µp0

√1− (x/a)2 + q′(x) (3.31)

Para satisfazer as condições na região de escorregamento (c < |x| < a), a q′ deve ser consideradaigual zero. Para achar essa variável na zona de adesão (|x| < c), considera-se o deslocamentorelativo igual zero (g′(x) = 0) e portanto obtém-se:

1

π=

∫ a

−a

q(ξ)

ξ − xdξ = 0 |x| ≤ c (3.32)

31

Substituindo a Equação (3.31) em (3.32) e resolvendo, tem-se:

q′(t) = −µp0c

a

√a− t2 |t| < 1 (3.33)

onde t = x/c . Para definir a meia largura da região de adesão, c, utiliza-se a equação de equilíbriotangencial da seguinte forma:

Q =

∫ a

−aq(ξ)dξ =

µp0π

2a(a2 − c2) (3.34)

De outra maneira, a Equação (3.34) pode ser escrita como:

c

a=√a− |Q/µP | (3.35)

As soluções apresentadas acima satisfazem as condições da zona de adesão e escorregamento.Essa solução também é conhecida como solução Mindlin-Cattaneo em homenagem a Cattaneo(1938) e Mindlin (1949) que publicaram seus trabalhos sobre efeito da força tangencial na interfacede contato de corpos elásticos. A Figura 3.8 mostra a influência da força de superfície tangencialna zona de adesão (Q/µP ).

As forças de superfícies explicadas anteriormente também podem ser apresentadas pela su-perposição de três curvas elípticas abaixo:

A distribuição de força normal máxima (−p0) em −a < x < a.

A distribuição de força tangencial máxima (µp0) em −a < x < a.

A distribuição de força tangencial máxima secundária (µp0c/a) em −c < x < c.

Considerando uma análise linear elástica, as superposições acima podem ser utilizadas paraobter as tensões e deformações do corpo. O código da solução Mindlin-Cattaneo implementado noprograma Matlab encontra-se no anexo A desse trabalho.

32

Figura 3.8: Tensões normalizadas de Mindlin-Cattaneo em regime de escorregamento parcial.

33

4 O Método dos Elementos de Contorno

4.1 A equação integral de contorno

A equação integral de contorno é obtida através da utilização do teorema de reciprocidadeda teoria da elasticidade, em conjunto com uma solução conhecida como solução fundamental daequação governante para uma carga discreta num corpo infinito.

O teorema de reciprocidade pode ser enunciado da seguinte forma:

"Se dois estados de equilíbrio distintos (u∗i , σ∗ij, b

∗i ) e (ui, σij, bi) existir em uma região de Ω+

na superfície limitada por Γ+, o trabalho realizado pelas forças de superfície e forças de corpo doprimeiro sistema sobre os deslocamentos do segundo é igual ao trabalho realizado pelas forças dosegundo sistema nos deslocamentos do primeiro."

O teorema pode ser provado usando o teorema da divergência e pode ser escrito como:

∫Γ

uiσ∗ijnjdΓ +

∫Ω

uib∗i dΩ =

∫Γ

u∗iσijnjdΓ +

∫Ω

u∗i bidΩ (4.1)

onde ui, σij e bi, são os deslocamentos, as tensões e as forças de corpo, e nj , são os componentesdo vetor normal, apontando para fora do domínio Ω. O domínio Ω∗ é denotado por um domínioinfinito e Ω é um domínio finito dentro de Ω∗ que tem um contorno Γ. As mesmas propriedadesdo material são assumidas por ambos os domínios. As forças de superfície Ti no contorno Γ sãodefinidas por:

Ti = σijnj (4.2)

Desta forma, a Equação (4.1) pode ser escrita como:∫Ω

(uib∗i − u∗i bi) dΩ =

∫Γ

(u∗iTi − uiT ∗i ) dΓ (4.3)

Os vetores de deslocamento u∗i , as forças de superfície T ∗i e as força do corpo b∗i na Equação (4.3)são escolhidos para ser a solução conhecida da equação de Navier devido a uma força pontual

35

unitária aplicada ao corpo, isto é:

µu∗i,kk +µ

1− 2νu∗k,ki + ∆(X− X′)ei = 0 (4.4)

onde ∆(X − X′) é a função delta de Dirac, XεΩ∗ é o ponto de carga singular e X′εΩ∗ é o pontocampo. O componente do vetor unitário ei em (4.4) corresponde a uma força positiva unitáriaaplicada em X′ na direção i. A função delta de Dirac tem a seguinte propriedade:∫

Ω∗g(X)∆(X− X′)dΩ(X) = g(X′) (4.5)

O campo de deslocamento e forças de superfície que correspondem à solução podem ser escritoscomo:

u∗i = U∗j δij = Uij(X′,X)ej (4.6)

eT ∗i = T ∗j δij = Tij(X′,X)ej (4.7)

onde δij é a função delta Kroneker. Uij e Tij são as soluções fundamentais dos deslocamentos eforças de superfícies, na direção j no ponto X devido a uma força pontual atuando na direção i emX′.

A componente de força de corpo b∗i (força por unidade de volume) corresponde a uma forçapontual e é dada por:

b∗i = ∆(X′,X)ei (4.8)

Substituindo b∗i na função de delta de Dirac, na Equação (4.3) e especificando X como a variávelde integração, obtém-se a seguinte equação considerando a força unitária atuando na direção i:∫

Ω

uj(X)∆ij(X′ − X)− Uij(X′,X)bj(X′)dΩ(X) =∫Γ

Uij(X′,X)Tj(X)− Tij(X′,X)uj(X)dΓ(X) (4.9)

Agora, usando as propriedades da função delta de Dirac, a Equação (4.9) resulta em:

36

ui(X′) =

∫Γ

Uij(X′ − x)Tj(x)dΓ(x)−∫

Γ

Tij(X′ − x)uj(x)dΓ(x)+∫Ω

Uij(X′,X)bj(X)dΩ(X) (4.10)

onde x′, x εΓ e X′,X εΩ. Esta é a famosa função Somigliana para deslocamentos, ui é uma repre-sentação contínua para os deslocamentos em qualquer ponto interior X ′ no domínio Ω. O campode tensão ao longo do corpo pode ser obtido através derivadas da Equação (4.10) como a seguir:

ui,k(X′) =

∫Γ

Uij,k(X′ − x)Tj(x)dΓ(x)−∫

Γ

Tij,k(X′ − x)uj(x)dΓ(x)+∫Ω

Uij,k(X′,X)bj(X)dΩ(X) (4.11)

A função Somigliana para tensões num ponto interior é obtida por substituição da Equação (4.11)na lei de Hooke (2.12), da qual obtém-se:

σij(X′) =

∫Γ

Uijk(X′ − x)Tk(x)dΓ(x)−∫

Γ

Tijk(X′ − x)uk(x)dΓ(x)+∫Ω

Uijk(X′,X)bk(X)dΩ(X) (4.12)

Uijk e Tijk são combinações lineares das derivadas das soluções fundamentais Uij e Tij eserão dadas na próxima seção.

4.2 Solução Fundamental

Para um problema no estado plano de deformação, as soluções fundamentais para desloca-mentos Uij(X′, x) e forças de superfície Tij(X′, x), definidas na Equação (4.5), são dadas por:

Uij(X′, x) =1

8πG(1− ν∗)(3− 4ν∗) ln(

1

R)δij +R,iR,j (4.13)

e

37

Tij(X′, x) =

−1

4π(1− ν∗)R[(1− 2ν∗)δij + 2R,iR,j]

∂R

∂n− (1− 2ν∗)(R,inj −R,jni) (4.14)

O campo de deformação fundamental Ukij(X′,X) e o campo de tensão Tkij(X′,X) comoapresentado na Equação (4.12) são dados por:

Ukij(X′,X) =1

4π(1− ν∗)R(1− 2ν∗)(R,jδki +R,iδkj +R,kδij) + 2R,iR,jR,k (4.15)

e

Tkij(X′,X) =µ

2π(1− ν∗)R22∂R∂n

[(1− 2ν∗)δijR,k + ν∗(R,jδik +R,iδjk)− 4R,iR,jR,k]

+ 2ν∗(niR,jR,k + njR,iR,k) + (1− 2ν∗)(2nkR,iR,j + njδik + niδjk)

− (1− 4ν∗)nkδij, i, j = 1, 2 (4.16)

Nestas equações, δij denota o delta de Kronecker, R(X′,X) representa a distância real entreo ponto fonte x′ e o ponto campo x, que é dada por:

R = |X− X′| e R,i =∂R

∂xi(4.17)

As soluções fundamentais para o estado de tensão plana podem ser obtida pela seguinte substituiçãoda relação modificada de Poisson e módulo de Young:

ν∗ =ν

1 + ν(4.18)

E∗ = E

[1− ν∗2

(1 + ν∗)2

](4.19)

38

4.3 A equação integral de contorno

A equação integral de contorno é obtida por um processo de limite, fazendo um ponto X′ nointerior do domínio Ω tender a um ponto x′ no contorno Γ. Este processo pode ser demonstradopela Figura 4.1 e a Equação (4.9) pode ser escrita como:

x

y

z

n

Ω

Γ

Γε'

Γε

ε

x'

Figura 4.1: O ponto fonte x′ pertence ao contorno.

ui(X′) =

∫Γ−Γe−Γ′e

Uij(X′, x)tj(x)dΓ(x)−

∫Γ−Γe−Γ′e

Tij(X′, x)uj(x)dΓ(x)+∫

Γ−Γe−Γ′e

Uij(X′, X)bj(X)dΩ(X) (4.20)

onde o contorno total é definido como:

Γ′ = (Γ− Γ′ε) + Γ′ε (4.21)

Γ′ε representa o contorno de um semicírculo de raio ε, Γ′ε tende a Γε quando ε → 0. Tomando olimite de ε→ 0, a formulação direta de elementos contornos é obtida como:

39

Cij(x′)uj(x

′) =

∫Γ

Uij(x′, x)tj(x)dΓ(x)−

∫Γ

Tij(x′, x)uj(x)dΓ(x)+∫

Ω

Uij(x′, X)bj(X)dΩ(X)dΩ(x) (4.22)

ondeCij(x

′) = δij + limε→0∫

Γε

Tij(x′, x)uj(x) (4.23)

A Equação (4.22) representa a formulação direta do método de elementos de contorno o qual re-laciona deslocamentos e forças de superfície no contorno. Esta equação integral de contorno paraum ponto geral sobre o contorno na ausência de forças do corpo bj pode ser escrita como:

Cij(x′)uj(x

′) +

∫Γ

Tij(x′, x)uj(x)dΓ(x) =

∫Γ

Uij(x′, x)tj(x)dΓ(x) (4.24)

onde Cij(X ′) = δij quando o ponto x′ está dentro do domínio Ω. A avaliação de Cij(x′) é maiscomplicada quando x′ está no contorno Γ, mas de um modo geral, tem-se:

Cij =

δij, se X′ ∈ ao domínioθint

2πδij, se X′ ∈ ao contorno

0, se X′ /∈ ao domínio ou ao contorno

(4.25)

onde θint é o ângulo interno do contorno Γ no ponto X′ (veja Figura 4.2).

Quando o ponto fonte encontra-se em ponto suave do contorno, isto é, não é um canto, tem-se:

Cij =1

2δij. (4.26)

Baseado na formulação apresentada, obtém-se a equação:

Cij(x′)uj(x

′) +

∫Γ

Tij(x′, x)uj(x)dΓ(x) =

∫Γ

Uij(x′, x)tj(x)dΓ(x) (4.27)

40

n

nX´

θint

Figura 4.2: Ângulo interno no contorno do componente.

Dividindo o contorno Γ em Ne elementos de contorno, a Equação (4.27), de outra forma,pode ser escrita como:

Cijuj(d) +Ne∑j=1

∫Γj

T ∗ikuidΓj =Ne∑j=1

∫Γj

U∗iktjdΓj (4.28)

Essa equação é aplicada em cada um dos nós do elemento de tal forma que a equação integralde contorno é transformada em um sistema linear de equações algébricas:

Hu = Gt (4.29)

onde as matrizes H e G contém as integrais das soluções fundamentais de forças de superfície Tije de deslocamentos Uij , e os vetores t e u contém todas as forças de superfícies e deslocamentosconhecidos ou não. Através de algumas manipulações algébricas podemos isolar as incógnitas emum vetor x de forma que o sistema (4.29) possa ser representado por:

Ax = b (4.30)

onde uma solução única pode ser obtida.

41

4.4 Elementos de contorno quadráticos contínuos

Na discretização utilizando elementos de contorno quadráticos contínuos, a geometria é apro-ximada por uma função quadrática ao longo de cada elemento, sendo necessários três pontos nodaispor elemento conforme mostrado na Figura (4.3).

Figura 4.3: Elementos quadráticos contínuos.

Assim, os deslocamento e as forças de superfícies são aproximados da seguinte forma:

u =

u1

u2

=

[N (1) 0 N (2) 0 N (3) 0

0 N (1) 0 N (2) 0 N (3)

]

u(1)1

u(1)2

u(2)1

u(2)2

u(3)1

u(3)2

= Nu(n) (4.31)

42

t =

t1

t2

=

[N (1) 0 N (2) 0 N (3) 0

0 N (1) 0 N (2) 0 N (3)

]

t(1)1

t(1)2

t(2)1

t(2)2

t(3)1

t(3)2

= Nt(n) (4.32)

onde u(n)i e t(n)

i são os valores nodais de deslocamentos e forças de superfícies, respectivamente, eN (i) são as funções de forma quadráticas definidas por:

N (1) =1

2ξ (ξ − 1) (4.33)

N (2) = 1− ξ2 (4.34)

N (3) =1

2ξ (ξ + 1) (4.35)

onde ξ representa uma coordenada adimensional ao longo do elemento, ilustrada na Figura 4.4.

Figura 4.4: Coordenada adimensional ξ para elementos quadráticos continuos.

43

Da mesma forma que os deslocamentos e forças de superfície, a geometria do problematambém é aproximada por elementos de contorno quadráticos contínuos:

Figura 4.5: Funções de forma de elementos quadráticos contínuos.

x1

x2

=

[N1 0 N2 0 N3 0

0 N1 0 N2 0 N3

]

x11

x12

x21

x22

x31

x32

(4.36)

Considere que o domínio tenha sido dividido em Ne elementos de contorno. Substituindo asEquações (4.31), (4.32) e (4.36) na Equação (4.28), tem-se:

clul +Ne∑j=1

∫Γj

T NdΓ

uj =

Ne∑j=1

∫Γj

U NdΓ

tj (4.37)

44

Chamando

∫Γj

UNdΓ = G (4.38)

e

∫Γj

TNdΓ = H (4.39)

tem-seN∑j=1

Hljuj =N∑j=1

Gljtj (4.40)

ou, na forma matricialHu = Gt (4.41)

4.4.1 Integração das matrizes H e G quando o ponto fonte não pertence ao ele-mento

A integração dos termos das matrizes H e G quando o ponto fonte não pertence aos elementosé uma integração regular que pode ser realizada usando, por exemplo, quadratura de Gauss. Estaintegração é descrita a seguir:

H(j) =

∫Γj

TlkN(j)dΓ =

∫ 1

−1

TlkN(j)|J |dξ (4.42)

G(j) =

∫Γj

UlkN(j)dΓ =

∫ 1

−1

UlkN(j)|J |dξ (4.43)

onde |J | representa o módulo do jacobiano da transformação (X1,X2)→ ξ:

45

|J | = dΓ

dξ=

(dx1

dξ

)2

+

(dx2

dξ

)21/2

(4.44)

onde dx1/dξ e dx2/dξ são obtidos derivando-se as equações (4.32) em relação a ξ.

Integrais singulares da ordem O(lnr) podem ser avaliadas eficientemente pela quadratura deGauss com uma transformação de variáveis cúbica que cancela exatamente a singularidade loga-rítmica. Uma outra possibilidade é o uso da quadratura logarítmica de Gauss. De acordo com estemétodo, os termos incluindo singularidades logarítmicas podem ser integrados por:

I =

∫ 1

0

ln

(1

ξ

)f(ξ)dξ ∼=

N∑i=1

wif(ξ) , (4.45)

onde N é o número de pontos de Gauss logarítmico.

Neste trabalho, os termos não singulares das matrizes H e G são integrados utilizando-sequadratura de Gauss padrão com 10 pontos de integração. Os termos singulares de G são do tipoln(r) sendo integrados usando quadratura logarítmica de Gauss com 10 pontos de integração. Já ostermos singulares de H são do tipo 1/r e precisam ser calculados no sentido do valor principal deCauchy. Uma maneira bastante simples de se tratar esta singularidade é através de considerações decorpos rígidos. Assumindo que um corpo rígido tenha todos os seus pontos do contorno deslocadosde um valor unitário e que não existam forças de corpo (bi = 0) na direção de um dos eixos decoordenadas, as forças de superfície em qualquer ponto do contorno deste corpo devem ser zero.Desta forma, a Equação (4.41) torna-se:

Hvq = 0 (4.46)

onde vq é um vetor que para todos os nós tem deslocamentos unitários ao longo da direção q e zerona outra direção. Para satisfazer a Equação (4.46) tem-se:

Hii = −Ne∑j=1

Hij j 6= i (4.47)

46

sendo j par ou ímpar.

O termo da diagonal da matriz H é igual a soma de todos os outros termos fora da diagonalcorrespondentes ao grau de liberdade em consideração.

4.5 Tensões no contorno

Para se calcular o tensor de tensões em um dado nó do contorno, considere um nó em que asdireções dos vetores tangente e normal ao contorno não coincidam com as direções dos eixos geo-métricos (Figura 4.6). Neste nó é criado um novo sistema de referência x′1x

′2 possuindo direções que

coincidam com os vetores tangente e normal ao contorno neste nó. Escrevendo os deslocamentos eas forças de superfícies neste sistema local tem-se:

u′i = lijuj

t′i = lijtj (4.48)

onde lij são os cossenos diretores.

No sistema local tem-se a seguinte relação:

σ′22 = t′2

σ′12 = t′1 (4.49)

A deformação ε′11 pode ser calculada, sabendo que:

47

σ'11

σ'12σ'22

Figura 4.6: Tensões no contorno.

ε′11 =1

2(u′1,1 + u′1,1) = u′1,1

u′1,1 =du′1dx′1

=du′1dξ

dξ

dx′1(4.50)

Usando geometria diferencial na Equação (4.50), pode-se notar que a direção local x′1 é tan-gente ao comprimento infinitesimal de arco ds dado por

ds =

√dx′1

2 + dx′22 =

√(dx′1dξ

)2

+

(dx′2dξ

)2

dξ

ds

dξ= J (4.51)

Um pequeno movimento ao longo de s corresponde a um pequeno movimento em x′1. Istopermite com que x′1 na Equação (4.50) seja substituído pela equação (4.51), ou seja:

48

ε′11 =du′1dξ

dξ

ds

ε′11 =du′1dξ

J−1 (4.52)

sendo

u1 =3∑i=1

N (i)u(i)1

du1

dξ=

3∑n=1

dN (i)

dξu

(i)1

(4.53)

onde N (i) são as funções de forma. Pode-se então obter a deformação

ε′11 =3∑

n=1

dN (i)

dξu

(i)1 J

−1 (4.54)

Na relação tensão deformação dada por:

σ′ij = λδijε′kk +Gε′ij (4.55)

tem-se três equações incógnitas σ′11, ε′22, ε

′12, que agora podem então ser calculadas.

Por último, as densidades de força tem que ser escritas no referencial global x1x2, ou seja

49

σ11

σ22

σ12

= T−1

σ′11

σ′22

σ′12

(4.56)

onde T é a matriz de transformação de coordenadas.

50

5 Método de elementos de contornos para problemas de contato

5.1 Introdução

O contorno de dois corpos isotrópicos homogêneos linearmente elásticos 1 e 2 são represen-tados por Γ1 e Γ2. Segundo Man (1994), quando dois corpos entram em contato, apenas uma partedo contorno de um corpo entra em contato com uma parte do contorno do outro. Portanto, os seuscontornos totais podem ser divididos em regiões de contorno em contato Γc e contorno sem contato

Γnc, tal como ilustrado na Figura 5.1.

Γ1 = Γ1nc + Γ1

c

Γ2 = Γ2nc + Γ2

c (5.1)

Figura 5.1: Região de contato Γc e região sem contato Γnc.

A região potencial de contato, Γc, pode conter regiões de adesão (Stick, st), regiões de escor-regamento (Slip, sl) e regiões de separação (sp). Assim, o contorno Γc para os corpos 1 e 2 pode serexpresso como:

Γ1,2c = Γ1,2

st + Γ1,2sl + Γ1,2

sp (5.2)

51

Supõe-se a força de atrito que atua ao longo da região de escorregamento segue a lei deatrito de Coulomb. Assim sendo, a solução numérica para os problemas de contato requer ummodelo matemático que possa representar os corpos em contato. Esta solução é representada pelaequação integral de contorno, formulada em termos das forças de superfície e deslocamentos nocontorno do domínio elástico. Em locais onde as forças de superfície são prescritas no contorno, osdeslocamentos correspondentes são as incógnitas da equação integral de contorno, e vice-versa.

Se dois corpos, sujeitos a uma carga externa, estão em contato sobre uma área Γc, a deforma-ção pode ser descrita por duas equações integrais acopladas, uma para cada corpo, como exposto aseguir:

C1iju

1j +

∫Γ1nc

T 1iju

1jdΓ1 +

∫Γ1c

T 1iju

1jdΓ1 =∫

Γ1nc

U1ijt

1jdΓ1 +

∫Γ1c

U1ijt

1jdΓ1 (5.3)

C2iju

2j +

∫Γ2nc

T 2iju

2jdΓ2 +

∫Γ2c

T 2iju

2jdΓ2 =∫

Γ2nc

U2ijt

2jdΓ2 +

∫Γ2c

U2ijt

2jdΓ2 (5.4)

Obtém-se uma solução numérica para o problema das equações integrais de contorno (5.3)e (5.4), discretizando separadamente os contornos dos corpos 1 e 2. Isto produz dois conjuntos deequações (um para 1 e outro para 2), dados por:

c1iju

1j +

N∑n=1

Hn1

ij un1

j =N∑n=1

Gn1

ij tn1

j (5.5)

c2iju

2j +

M∑m=1

HmB

ij um2

j =M∑m=1

Gm2

ij tm2

j (5.6)

onde N e M são o número total de nós dos corpos 1 e 2 respectivamente. Dessa forma, doisconjuntos de equações lineares são obtidos. Este conjunto pode ser expresso na forma de uma

52

matricial como:H1u2 = G1t1 e H2u2 = G2t2 (5.7)

Os vetores de u1,2 e t1,2 representam, respectivamente, valores de deslocamentos e forçasde superfície nos contornos dos corpos 1 e 2. Na região de contato, os dois sistemas de equaçõescompartilham as variáveis de contorno do problema, ou seja, as equações são acopladas e devemser resolvidas simultaneamente para qualquer combinação de carga externa e qualquer condição decontato. Se as condições de contorno, dentro e fora de qualquer região de contato, são implemen-tadas, as equações (5.7) podem ser reduzidas a um sistema unificado:

Ax = b (5.8)

Para problemas lineares, uma vez que este sistema de equações foi resolvido, pode ser ob-tida uma solução final para os deslocamentos e forças de superfície em todos os lugares sobre oscontornos. Entretanto, os problemas de contato podem ser não lineares e a extensão da zona decontato pode ser desconhecida. Neste caso, a determinação da zona de contato é parte da soluçãodo problema. Isto significa que alguns problemas de contato exigem um procedimento iterativo desolução. Durante o processo iterativo, coeficientes de A e b, derivados a partir do interior da zonade contato, podem ser alterados de uma iteração para outra. O número de mudanças na matriz Aé pequeno, dado que o número de elementos que necessitam alterações nas condições de contatoé geralmente uma fração pequena do total. Neste caso, o sistema de equações tem de ser reorde-nado para a próxima iteração, de modo a acomodar as alterações na zona de contato para que, emseguida, a matriz A seja atualizada. A repetição deste procedimento, até que a solução final sejaencontrada, é ineficiente e cara.

Abascal (1995) sugere resolver o sistema de equações atualizado de forma eficiente e semrecorrer a uma reformulação da matriz de todo o sistema. Ele sugere as variáveis desconhecidasna zona de contato potencial separadas das incógnitas fora dela. Esta técnica pode acelerar consi-deravelmente a solução iterativa. No entanto, deve-se considerar que uma vez que a atual zona decontato pode ser desconhecida, é essencial que uma zona de potencial contato a ser escolhida sejamaior do que a região de contato final. Para as zonas de contato potenciais, as equações obtidas apartir das condições de contato tem que ser expressas de forma explícita, de modo que elas possamser separadas daquelas que estão fora da zona de contato.

53

Se os números totais de nós fora da zona de contato potencial dos corpos 1 e 2 são N1nc e

M2nc respectivamente, isto resulta em 2(N1

nc + M2nc) equações lineares, uma vez que existem duas

incógnitas por nó. No interior da zona de contato em cada nó de contorno, ambas as componentes deforça de superfície e os dois componentes de deslocamento são desconhecidos. Assim, para um parde nós de contato existem oito incógnitas. Para calcular estas oito incógnitas, quatro equações sãoencontradas ao considerar a compatibilidade de deslocamento e equilíbrio de forças de superfíciena interface de contato. As outras quatro equações são obtidas a partir de equações integrais decontorno. Essas equações de compatibilidade e de equilíbrio podem ser obtidas explicitamentepara cada par de nós com potencial de contato, considerando o estado do contato do par do nó emsi com seus pares de nós vizinhos imediatos.

A partir disso, um sistema de equações de contato pode ser obtido na forma de matriz, demodo que ele possa ser facilmente incorporado ao conjunto da matriz final. Além disso é necessárioconsiderar que a região de contato pode conter uma combinação de separação, adesão e desliza-mento de pares de nós, sendo que a ordem na qual eles ocorrem depende do problema envolvido.O sistema de equações de contato pode ser formulado para lidar com qualquer situação possível.

Este sistema de equações é escrito em termos de um sistema de coordenadas locais com asnormais unitárias exteriores tomadas como sendo positivas. Para todos os pares de elementos decontato potenciais, as equações podem ser descritas como condições de contorno dentro da regiãode contato. Portanto, a matriz A tem uma parte linear onde consta as equações integrais de contorno(AL), e outra parte não linear que contém as equações de contato (ANL). As equações de contato sãonão lineares, pois podem mudar de um incremento para outro. Assim sendo, as duas partes devemser resolvidos simultaneamente, pois todas as equações da parte linear são indefinidas e precisamde dados das equações não lineares para serem resolvidas. Dessa forma, a matriz A tem a seguinteforma:

A =

[AL

ANL

](5.9)

onde:

54

AL =

[A1nc 0 H1

c 0 -G1c 0

0 A2nc 0 H2

c 0 -G2c

](5.10)

e

ANL =[

0 0 C1u C2

u C1t C2

t

](5.11)

A matriz AL é calculada apenas uma vez pois seus elementos não variam com as condiçõesde contato nem com o incremento do carregamento (desde que as deformações se mantenhampequenas). Incorporando as matrizes (5.10) e (5.11) na matriz principal, tem-se:

A1nc 0 H1

c 0 -G1c 0

0 A2nc 0 H2

c 0 -G2c

0 0 C1u C2

u C1t C2

t

x1nc

x2nc

u1c

u2c

t1c

t2c

=

b1

b2

v1,2

(5.12)

As matrizes A1nc e A2

nc são as matrizes A dos corpos 1 e 2, respectivamente, considerandoapenas os nós que não estão na região de potencial contato; as matrizes H1

c e H2c são as matrizes H

dos corpos 1 e 2, respectivamente, considerando os nós na região de potencial contato; as matrizesG1c e G2

c são as matrizes G para os corpos 1 e 2, respectivamente para os nós na região de poten-cial contato; as matrizes C1

u e C2u representam as matrizes de compatibilidade de deslocamentos

considerando as condições de contato dos nós dos corpos 1 e 2, respectivamente; os vetores x1nc

e x2nc contém os graus de linearidade desconhecidos na região de não contato dos corpos 1 e 2,

respectivamente; u1c e u2

c são os vetores de deslocamentos dos corpos 1 e 2 da região de contato,respectivamente; t1

c e t2c são os vetores de forças de superfície dos corpos 1 e 2, respectivamente, na

região de contato; b2 e b2 são os vetores b para os corpos 1 e 2, respectivamente e, finalmente, v1,2

é um vetor que contém a distância entre dois nós em contato (g0).

As matrizes Cu e Ct são construídas com as matrizes cu e ct para cada par de nó na região

55

possível de contato. Seus valores dependem da situação de contato entre os nós e são dados por:

Escorregamento:

cu =

0 0 0 0

0 0 0 0

1 0 1 0

0 0 0 0

ct =

1 0 −1 0

±µ 1 0 0

0 0 0 0

0 0 ±µ 1

(5.13)

Adesão:

cu =

0 0 0 0

0 0 0 0

1 0 1 0

0 1 0 1

ct =

1 0 −1 0

0 1 0 −1

0 0 0 0

0 0 0 0

(5.14)

Separação:

cu =

0 0 0 0

0 0 0 0

0 0 0 0

0 0 0 0

ct =

1 0 0 0

0 1 0 0

0 0 1 0

0 0 0 0

(5.15)

Para se resolver o sistema não linear da Equação (5.12), Rodríguez-Tembleque Solano (2009)sugeriu o uso de uma função vetorial residual, dada por:

R = Ax− b (5.16)

O gradiente analítico de R é conhecido e é dado por:

56

∇R = A (5.17)

Deve-se então encontrar o vetor x para o qual R = 0. Existem vários métodos numéricosusados para encontrar o zero de função não linear. Dentre eles, o método de Newton, está entreos mais usados e é o método usado neste trabalho. O método de solução implementado mostrouboa eficiência e convergiu em menos de 12 iterações em cada incremento nas condições específicasconsideradas neste trabalho.

5.2 Aplicação incremental de carga

Nos problemas de contato, a solução final é obtida somente quando a carga externa máximaou contato total da região de contato é atingida. Entretanto, como a região potencial de contatoé uma função da força aplicada e do coeficiente de atrito, a solução deve ser obtida com auxíliode uma solução numérica iterativa. Ainda mais se o contato for entre contornos não conformes,a solução aproximada é alcançada com aplicação gradual de carga em várias etapas. É necessárioconsiderar também que nos problemas de contato, a equação de equilíbrio e compatibilidade aolongo da região de contato devem ser satisfeitos.

Sabe-se que existem vários métodos para a aplicação de carga nos problemas de contato.Aplicação instantânea de uma carga externa em uma única vez é aceitável sé houver contato nó anó na região de adesão. Nesse caso, o processo iterativo de solução do sistema é utilizado até que aregião de escorregamento e adesão sejam conhecidos.

Outra abordagem, mais recomendada, é a aplicação incremental de carga. Nesse caso, a cargatotal será aplicada gradualmente, em vários incrementos. A vantagem dessa abordagem é um resul-tado mais aproximado e também sua aplicação, uma vez que ela serve tanto para modo de contatonó a nó quanto para contato de uma superfície com um nó. Por ser um processo iterativo, a históriadas tensões de contato, as quais tem importância em casos específicos de contato como fretting,devem ser registradas.

Na abordagem incremental, as condições de contorno, geometria deformada dos corpos eestado de contato são atualizados em cada iteração. O número de incrementos é predefinido pelo

57

usuário e a soma dos incrementos é igual à carga total externa, ou seja:

Pmax =M∑m=1

∆Pm (5.18)

onde ∆P é a carga aplicada em cada incremento, Pmax é carga externa total aplicada eM é númerode incrementos. A força total aplicada é calculado pela Equação (5.19).

Pmj = Pm−1

j + ∆Pmj (5.19)

onde Pmj é a força externa total até o incremento m. Na aplicação de cada incremento de carga

(∆P ) haverá uma resposta do sistema que se encontra inicialmente em equilíbrio. A resposta dosistema resulta em pequenas variações nas forças de superfícies e nos deslocamentos no contornodo sistema até obter um novo equilíbrio. As variações de deslocamento e forças de superfície sãodefinidas como:

umj = um−1j + ∆umj

tmj = tm−1j + ∆tmj (5.20)

onde ∆umj e ∆tmj são variações incrementais de deslocamento e força de superfície devido a cargaincremental ∆Pm

j . Rescrevendo a equação integral de contorno, tem-se:

Cij(um−1j + ∆umj ) +

∫Γ

Tij(um−1j + ∆umj )dΓ =

∫Γ

Uij(tm−1j + ∆tmj )dΓ

A integral do contorno para cada incremento é obtida através da equação abaixo:

Cij∆umj +

∫Γ

Tij∆umj dΓ =

∫Γ

Uij∆tmj dΓ

58

Os nós do contorno discretizado devem ser organizados na região de contato e formar paresde nós, criando a superfície de contato. Todos os pares de nós devem ser resolvidos como um sis-tema de contato independente, satisfazendo forças de superfície incremental ∆t e deslocamentosincrementais ∆u. Além disso este procedimento deve ser aplicado tanto para a equação de equi-líbrio quanto para a equação de compatibilidade. A partir disso, percebe-se que para cada corpohaverá um sistema de equações, mostrado anteriormente na Equação (5.7). Reformulando para ocaso incremental, tem-se:

H1∆u1 = G1∆t1 e H2∆u2 = G2∆t2 (5.21)

onde os vetores ∆u e ∆t contém, respectivamente, os valores de deslocamentos e forças de su-perfície para cada incremento. Para obter a solução do problema as condições de contorno fora daregião de contato e restrições de contato dentro da região de contato deverão ser aplicadas.

5.3 Sistema de coordenadas na região de contato

Nos problemas de contato, as variáveis fora da região de contato são representadas nas co-ordenadas globais, definidas pelas coordenadas cartesianas, enquanto as variáveis dentro da zonade contato são analisadas nas coordenadas locais. Portanto, será necessária uma transformação decoordenadas globais para coordenadas locais, na região potencial de contato, para todos os paresde nós. Essa transformação é feita através de rotação na direção média do vetor normal.

Considera-se os corpos A e B tendo como pares de nós a e b com coordenadas (xa, ya)

e (xb, yb) na zona de contato com vetor normal n e vetor tangente s, respectivamente, da formaapresentada na Figura 5.2.

Para encontrar a direção média do vetor normal utiliza-se a seguinte fórmula:

na =na − nb

|na − nb|= −nb

59

A

B

P

x

y

Figura 5.2: Sistema de coordenadas globais.

A distância entre os nós é determinada pela distância normal inicial entre os corpos, mostradona Figura 5.3.

A distância normal entre os nós, chamada de g0, é obtida pela Equação (5.22):

g0 =∣∣(yb − yb)∣∣ .cos(θ) +

∣∣(xb − xb)∣∣ .sen(θ)

= |dy| .cos(θ) + |dx| .sen(θ)

= |dy| nay.+ |dx| .nax (5.22)

onde nay e nax são componentes normal e tangencial da média do vetor normal, considerada anterior-mente. Sob carregamento externo, os corpos se aproximam e suas superfícies começam a ter contatoà medida que g0 vai diminuindo. Fisicamente tem-se g0 ≥ 0 pois a penetração dos corpos não serápermitida e o contato somente acontece quando g0 = 0. Por outro lado, matematicamente obter a

60

a

b

ña

s a

ñb

s b

sa

na

nbsb

ñnb

nag0

Figura 5.3: Sistema de coordenadas locais.

distância exatamente igual a zero é impossível, e portanto no programa implementado deve-se con-siderar uma faixa de valores pequenos, próximos à zero, para poder calcular a região de contato, ouseja, definir um zero flutuante. Vale reforçar que as forças de superfícies normais dos corpos (tan etbn) na região de contato devem obter valores negativos pois estão em compressão, ou seja:

tan = tbn ≤ 0

5.4 Restrições na região de contato

Na abordagem incremental, as condições de contato para os pares de nós dos corpos sãoexpressas relativas à carga incremental. Nesse sentido, ∆u e ∆t, obtidos em cada incremento,devem estar em equilíbrio de forças e estarem compatíveis com deslocamento total para cada nó.As equações de equilíbrio de força, de superfície e compatibilidade de deslocamento para cadamodo de contato, são listadas abaixo. Nessas expressões, o lado direito das equações representa t eu total, aplicadas até o incremento anterior, sendo portanto valores já conhecidos. No lado esquerdodas equações estão os valores de variações de t e u a serem obtidos nas coordenadas locais.

61

5.4.1 Adesão ou escorregamento parcial

∆(tat )m −∆(tbt)

m = −[(tat )

m−1 − (tbt)m−1]

∆(tan)m −∆(tbn)m = −[(tan)m−1 − (tbn)m−1

]∆(uat )

m −∆(ubt)m = 0

∆(uan)m −∆(ubn)m = g0 −[(uan)m−1 − (ubn)m−1

]= gm0 (5.23)

5.4.2 Escorregamento

∆(tat )m −∆(tbt)

m = −[(tat )

m−1 − (tbt)m−1]

∆(tan)m −∆(tbn)m = −[(tan)m−1 − (tbn)m−1

]∆(tat )

m − µ∆(tbt)m = −

[(tat )

m−1 ± µ(tbt)m−1]

∆(uan)m −∆(ubn)m = g0 −[(uan)m−1 − (ubn)m−1

]= gm0 (5.24)

5.4.3 Separação

∆(tat )m −∆(tbt)

m = −[(tat )

m−1 − (tbt)m−1]

∆(tan)m −∆(tbn)m = −[(tan)m−1 − (tbn)m−1

]∆(tat )

m = −[(tat )

m−1]

∆(tan)m = −[(tan)m−1

](5.25)

62

5.5 Decisão do estado de contato

A verificação do estado de contato, em qualquer fase, é baseada nas decisões de contatomostradas a seguir para um par de nós de contato a e b. Primeiro deve-se verificar se os paresseparados permanecem separados ou se entram em contato e vice-versa. Deve-se notar tambémque a violação representa uma incompatibilidade geométrica e não deve ocorrer em qualquer fasedos cálculos.

Nas fórmulas a seguir, t e u representam força de superfície e deslocamento e ∆t e ∆u sãoas variações incrementais deles. Os subscritos n e t são referentes à direção normal e tangencial,m− 1 representa as condições no passo anterior e m denota a situação depois da aplicação de umacarga adicional.

5.5.1 Separação ou Contato

Considera-se que os corpos separados continuam separados, se ocorrer:

(∆uan + ∆ubn)m < gm−10

sendo g0 a distância entre nós correspondentes dos corpos A e B. Caso contrário, os corpos entramem contato, ou seja:

(∆uan + ∆ubn)m > gm−10

De outra forma, considera-se que os corpos em contato continuam em contato caso:

tm−1n + ∆tmn < 0

63

caso contrário, os corpos perdem o contato, ou seja:

tm−1n + ∆tmn > 0

5.5.2 Adesão ou Deslizamento

Depois de conhecer o estado de separação ou contato entre os corpos, se ocorrer contato,deve-se aplicar as fórmulas a seguir para diferenciar a adesão e o deslizamento para aqueles paresde nós que estão em contato. A detecção é feita por violação de equilíbrio de força tangencial enormal. Essa violação geralmente implica que a continuidade de força sobre a região de contatofoi violada. Neste caso, uma redistribuição da força é necessária. Faz-se necessário lembrar quequalquer redistribuição de força no interior da região de contato é alcançada por meio da variaçãodas dimensões das regiões de adesão e deslizamento.

Nesse caso, basta considerar os corpos no estado de adesão e verificar a validade da fórmulaabaixo:

|tmt | = |tm−1t + ∆tmt | 6 µ|tm−1

n + ∆tmn | = µ |tmn | (5.26)

Se a condição 5.26 for válida, os corpos continuam no estado de adesão. Caso contrário, se tiver:

|tm−1t + ∆tmt | > µ|tm−1

n + ∆tmn |

os corpos vão deslizar um sobre o outro.

Pela lei de atrito de Coulomb, para o coeficiente de atrito entre dois corpos existe um co-eficiente de atrito estático (µs), e um coeficiente dinâmico (µd). Apesar que a diferença entre oscoeficientes é muito pequena, a consideração de somente um coeficiente de atrito seria idealizaçãodo fenômeno. Todavia, nos problemas onde ocorre regiões de adesão e escorregamento simulta-neamente, o µ é aproximado para um valor médio (Hills e Nowell, 1994). Nos testes feitos nolaboratório do Oxford, esse valor médio para ligas de alumínio sob condições de fadiga por fretting

é baixo, da ordem de 0,2.

64

6 Método de solução de sistemas não lineares

Sistemas não lineares são frequentemente encontrados nos problemas reais de física e mate-mática. Pelo aspecto físico, sistemas não lineares são sistemas que não satisfazem o princípio desuperposição, ou seja, a saída do sistema não é linearmente proporcional à entrada do sistema. Poroutro lado, para matemáticos, sistemas não lineares são sistemas que o efeito de fatores externosnão é puramente aditivo e não alteram proporcionalmente a uma mudança de uma entrada.

Esse tipo de sistema geralmente é aproximado por sistemas lineares e, dependendo do mé-todo aplicado e modo de aplicação, obtém-se uma solução próxima da realidade. Entretanto, ossistemas não lineares não podem ser escritos em forma de combinação de sistemas lineares exatos.Além disso, essa aproximação em alguns casos não consegue demonstrar o real comportamento dosistema, como por exemplo, na singularidade ou teoria de caos.

Existem vários métodos para resolver sistemas não lineares como secante, quase-Newton,Newton, entre outros. Aqui somente o método de Newton será abordado, pois este foi o métodoutilizado para resolver o sistema não linear de contato.

6.1 Princípios do método de Newton

O método de Newton (também conhecido como Newton-Raphson) usa a derivada da função,repetindo em várias iterações até achar o ponto mínimo da função (Figura 6.1). O uso desse métodoé mais comum nas funções onde a obtenção de derivada é acessível. O objetivo do método deNewton é encontrar o valor x∗ no qual se tem f(x∗) = 0. Aproximando a função desejada por umasérie de Taylor em torno de xk, tem-se:

f(x∗) ' f(xk) + (x∗ − xk)f ′(xk) (6.1)

Para se achar o zero da função, tem-se a fórmula seguinte:

65

Figura 6.1: Ilustração da iteração do método de Newton.

f(xk) + (x∗ − xk)f ′(xk) = 0 (6.2)

A Equação (6.2) pode ser escrita da seguinte forma:

xk+1 = xk −f(xk)

f ′(xk)= xk + δk (6.3)

onde xk+1 seria uma nova estimativa para o sistema. A formulação dada precisa de um ponto inicialpara x0 e dois valores pequenos maiores que zero (ε1, ε2 > 0) como critério para interromper oprocesso iterativo. Os critérios são definidos da forma seguinte:

|xk+1 − xk| < ε1(1 + |xk+1|)

|f(xk+1)| < ε2 (6.4)

66

Se a derivada da função não for conhecida, os métodos numéricos podem auxiliar para teruma derivada numérica da função. Por outro lado, no domínio considerado, a derivada da funçãonão pode ser igual zero, se não o método fica preso num ponto de mínimo local. Além disso, comoestá ilustrado na Figura 6.2, o método pode ficar preso num ciclo entre dois pontos com derivadasiguais.

Figura 6.2: Demostração de caso de comportamento cíclico.

Para resolver esse problema, é sugerido usar o método de Newton amortecido (Buchanan eFitzgibbon, 2005), qual substitui a iteração normal por:

xk+1 = xk +δk2j

(6.5)

onde

j = min i : 0 ≤ j ≤ imax∣∣∣∣f (xk − δk

2jf(xk)

f ′(xk)

)∣∣∣∣ < |f(xk)| (6.6)

67

Caso esta condição seja impossível de se cumprir, pode-se considerar j = 0. Usualmente,utiliza-se imax = 4, mas esse valor pode ser ajustado dependendo do problema. Valores maiorespara imax pode aumentar o custo computacional.

6.2 Sistema multidimensional

Suponha um sistema não linear de n equações. Isso significa calcular n variáveis simultane-amente.

fn(x1, ....., xn) = 0 (6.7)

Considerando o caso unidimensional na seção anterior e o ampliando para caso multidimen-sional, tem-se:

F (xk) +∇F (xk)(x∗ − xk) = 0 (6.8)

E o ponto para cada iteração seria:

xk+1 = xk + δk (6.9)

onde

δk = (∇F (xk))−1F (xk) (6.10)

As mesmas considerações feitas para o caso unidimensional podem ser aplicados aqui. Paramétodo de Newton amortecido, tem-se a equação seguinte:

68

j = min i : 0 ≤ j ≤ imax∥∥∥∥f (xk − δk

2jf(xk)

f ′(xk)

)∥∥∥∥2

< ‖f(xk)‖2 (6.11)

O critério de parada utilizado nesse trabalho é a variação do valor de x para ε1 = 10−8, ouseja:

|xk+1 − xk| < 10−8(1 + |xk+1|) (6.12)

O método de Newton implementado nesse trabalho não convergiu para valores elevados deforça tangencial. O autor acredita que isso é devido de deslocamentos maiores entre os corpos ecausa separação completo dos corpos.

6.3 Algoritmo do programa principal

1. Entrada de dados

2. Para i = 1, 2 (monta as matrizes para os corpos 1 e 2)

(a) Calcule as matrizes Gi e Hi.

(b) Aplicar as condições de contorno na região de não contato (gera as matrizes Ai e bi).

(c) Aloca a matriz Ai e o vetor Hi na matriz global ANL e bNL.

3. u0 = 0

4. t0 = 0

5. Para i = 1 até número de passos de carga:

(a) Calcula a distância entre os nós na região de provável contato.

(b) Aplica o método de Newton para resolver o sistema não linear R = Ax− b = 0.

(c) Monta os vetores ∆u e ∆t baseado nas condições de contorno e nas variáveis desco-nhecidas x.

(d) ui = ui−1 + ∆u

69

(e) ti = ti−1 + ∆t

6. Fim do programa.

6.4 Algoritmo para método de Newton

Repita:

1. Monte as matrizes ANL e bNL.

2. Monte as matrizes A e b.

3. R = Ax− b

4. d = −A−1.R

5. x = x0 + d

6. δ =√

(x− x0)′.(x− x0)

7. x0 = x

Até que δ < ε

6.5 Algoritmo para a montagem das matrizes ANL e bNL

Para i = 1 até número de pares de nos na região de possível contato

· Se |u1n + u2

n − h| < ε ou tn < 0 os nós não estão em contato.Monte cuij e ctij considerando os nós livre de contato.

· Se não, se |tt| ≥ tn, os nós estão em contato com deslizamento.Monte cuij e ctij considerando nós em contato deslizante.

· Se não, os nós estão em condição de adesão.Monte cuij e ctij considerando o contato em adesão.

70

7 Simulação numérica e resultados

Neste capítulo serão abordados os resultados obtidos com o programa implementado paraproblemas clássicos de contato. O problema geral analisado é o contato de dois cilindros com eixosparalelos ilustrada na Figura (7.1). Estes são validados com resultados analíticos encontrados naliteratura. Em cada exemplo, são apresentados dois gráficos com número de elementos diferentesde maneira que será possível comparar a influência do refinamento da malha utilizada. Em todosos exemplos foram utilizados elementos de contorno quadráticos contínuos.

d

Figura 7.1: Contato de dois cilindros sob carregamento vertical e horizontal.

7.1 Exemplo 1: Contato de duas sapatas elásticas sem atrito sob carregamentovertical

O problema de contato entre dois cilindros foi aproximado pelo contato de duas sapatas, tendoem vista diminuir o custo computacional e simplificar o problema sem afetar o campo de interesse(Figura 7.2). Foram consideradas duas sapatas elásticas idênticas sob carregamento vertical P .Os dois corpos possuem diâmetro de 70 mm, módulo de elasticidade de 73,4 GPa, coeficientede Poisson de 0,33 e carga vertical de P =100 N/mm. Na região de contato, o campo de tensãodeste problema é bastante próximo do que ocorre no contato entre dois cilindros (Figura 7.2).

71

As condições de contorno consideradas são as mesmas utilizadas na literatura para problemas defretting, como por exemplo no trabalho do McVeigh et al. (1999).

Corpo A

Corpo B

p

2h

2w

Figura 7.2: Modelo de duas sapatas em contato.

Como explicado anteriormente no capítulo 5, quando dois corpos entram em contato, apenasuma parte do contorno de cada corpo entra em contato efetivamente. O modelo implementado dassapatas é ilustrado na Figura 7.3. Foi assumido w=h=6,5 mm.

Conforme a Figura 7.3, a sapata de baixo (sapata 2) está fixada nas direções y e fixada nadireção x no nó do meio da aresta para evitar a rotação da sapata e a sapata de cima (sapata 1) estásob carregamento vertical P . Este problema de contato de dois cilindros sob carregamento verticalsem atrito é um caso bastante clássico da literatura e é amplamente utilizado para validar o modelonumérico, pois possui solução analítica. Os resultados obtidos são comparados e validados a partirdesta solução.

72

−10 −5 0 5 10−8

−6

−4

−2

0

2

4

6

8

x (mm)

y (m

m)

Figura 7.3: Modelagem do problema no programa implementado.

7.1.1 Resultados numéricos

Foram utilizadas duas malhas, com números de elementos diferentes, para verificar a in-fluência delas na solução do problema. A Figura 7.4 mostra a tensão normal na região de contatoutilizando 60 elementos na aresta de contato, 12 elementos na aresta superior e 6 elementos emcada aresta lateral. A partir disso, a solução numérica é comparada com a solução analítica deHertz. Considera-se ainda que na ausência de atrito, a tensão de cisalhamento é igual a zero e porisso não foi apresentada na figura.

A Figura 7.4, com a malha menos refinada, já mostra uma boa concordância entre os resul-tados numéricos e a solução analítica. Todavia, a Figura 7.5, obtida com o dobro de elementos emcada face do corpo, mostra um resultado mais aproximado, especificamente no começo da regiãode contato.

Na Figura 7.6 se encontram duas curvas para cada corpo que representam deslocamentosao longo dos eixos x e y. O deslocamento ao longo do eixo x para o corpo 1 é positivo quandoos valores de x são negativos, zero no meio da região de contato e negativo para valores de x

73

−8 −6 −4 −2 0 2 4 6 80

100

200

300

400

500

600

700

x(mm)

Ten

sões

(M

Pa)

σ

y

σy analítico

Figura 7.4: Tensão normal com 60 elementos quadráticos contínuos na região de contato compa-rando com a solução de Hertz no contato sem atrito.

Figura 7.5: Tensão normal com com 120 elementos quadráticos contínuos na região de contatocomparando com solução de Hertz no contato sem atrito.

74

positivos. Isso significa que o corpo 1 sob carregamento vertical comprime a aresta de contato. Umcomportamento parecido, porém com menor magnitude, é obtido para o corpo 2.

−8 −6 −4 −2 0 2 4 6 8−0.06

−0.05

−0.04

−0.03

−0.02

−0.01

0

0.01

x (mm)

Des

loca

men

tos

(mm

)

uy corpo 1

ux corpo 1

uy corpo 2

ux corpo 2

Figura 7.6: Deslocamentos dos corpos ao longo do eixo x e y no contato sem atrito.

A Figura 7.6 mostra também os deslocamentos ao longo do eixo y. Pode-se ver que as extre-midades da região de contato do corpo 1 se desloca mais que a região central desse corpo. No corpo2, o deslocamento em y tem comportamento contrário, ou seja, se desloca mais no centro da regiãode contato que nas extremidades. No entanto percebe-se que o deslocamento dos corpos convergempara um valor igual na região central dos corpos. Não foram encontrados gráficos de deslocamentona literatura, mas os resultados apresentados são relevantes para quem estuda problemas de contato.

A Figura 7.7 mostra as tensões equivalentes de von Mises nos corpos 1 e 2. Note-se que atensão máxima na região de contato ocorre no interior dos corpos. Este resultado é bem conhecidodos livros de projetos de máquinas (Shigley, 1972) sendo que este fato é responsável pela nucleaçãoda trinca no interior de componentes mecânicos como pista de rolamentos e dentes de engrenagens.Esta trinca é a causa da falha por fadiga sub-superficial.

Para verificar as tensões nos pontos internos dos corpos ao longo da linha x = 0 existeuma solução analítica de Hertz, apresentada na seção 3.4. A Figura 7.8 mostra a comparação das

75

Figura 7.7: Mapa de cor das tensões equivalente de von Mises nos corpos sob carregamento verticalsem atrito.

tensões nos pontos internos com a solução analítica existente. A partir destes dados destaca-se queo resultado está em boa concordância com a solução analítica. O autor acredita que este erro podeser causado pelo uso de derivada das funções de forma, ao invés do uso de equação integral, paracalcular as tensões.

76

−700 −600 −500 −400 −300 −200 −100 00

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Tensão (MPa)

y (m

m)

σ

y analítico

σx analítico

τ analítico

σy numérico

σx numérico

τ numérico

Figura 7.8: Variação de tensões principais no plano ao longo do eixo y e comparação com a soluçãoanalítica.

77

7.2 Exemplo 2: Contato de duas sapatas elásticas com atrito sob carregamentovertical

Neste exemplo, o mesmo problema anterior foi considerado, com as mesmas dimensões epropriedades dos materiais, mas com atrito entre os corpos, como acontece nos casos reais. Naimplementação foi considerado um coeficiente de atrito µ = 0.2 usando a lei de atrito de Coulomb.

7.2.1 Resultados numéricos

A diferença principal entre o problema sem atrito e o com atrito é a tensão de cisalhamento.No caso do contato sem atrito, não existe tensão de cisalhamento, pois não existem forças tangenci-ais. No caso de contato com atrito, o próprio contato gera forças tangenciais e, consequentemente,tem a tensão de cisalhamento. As tensões obtidas com a malha mais refinada são apresentadas naFigura 7.9. A malha utilizada é a mesma do exemplo anterior.

−8 −6 −4 −2 0 2 4 6 8−100

0

100

200

300

400

500

600

700

x(mm)

Ten

sões

(M

Pa)

σ

y

τxy

/µ

σy analítico

Figura 7.9: Tensões normal e cisalhante no corpo 1 com 60 elementos quadráticos contínuos naregião de contato considerando µ = 0, 2.

Devido à simetria do problema, a tensão de cisalhamento é nula no ponto central de contato

78

(x = 0) e aumenta no sentido das extremidades da região de contato. Até certo ponto quando acondição de contato muda de adesão para deslizamento. Neste ponto ocorre um pico onde a tensãode cisalhamento passe a decrescer até a extremidade da região de contato onde é zero. Este resultadofoi obtido em apenas um incremento, ou seja, toda a carga vertical foi aplicada em uma só passo.Observa-se na Figura 7.9 que ao aplicar o carregamento vertical gradualmente (nesse caso em 50incrementos), a variação da tensão de cisalhamento na região de contato tem uma forma suave (sempicos).

−1.5 −1 −0.5 0 0.5 1 1.50

100

200

300

400

500

600

700

x (mm)

Ten

sões

(M

Pa)

ty (50 passos)

tx/µ (50 passos)

ty (1 passo)

tx/µ (1 passo)

Figura 7.10: Tensões normal e cisalhante no corpo 1 com a malha mais refinada no contato comatrito .

Para esse caso não existe solução analítica e não foi encontrado nenhum resultado numéricona literatura para validar a curva de τxy. O deslocamento dos corpos na direção y é igual do des-locamento obtido para caso sem atrito, entretanto o deslocamento na direção x altera um pouco.Como se pode observar na Figura 7.11 o deslocamento tangencial dos corpos ao redor do pontocentral de contato são idênticos, significando a adesão dos corpos.

O mapa de cor das tensões equivalentes de von Mises nos corpos para o caso de contato comatrito é ilustrado na Figura 7.12.

Observe-se pouca diferença entre os casos com e sem atrito para as tensões de von Mises.

79

−8 −6 −4 −2 0 2 4 6 8−0.06

−0.05

−0.04

−0.03

−0.02

−0.01

0

0.01

x (mm)

Des

loca

men

tos

(mm

)

uy corpo 1

ux corpo 1

uy corpo 2

ux corpo 2

Figura 7.11: Deslocamentos dos corpos ao longo do eixo x e y no contato com atrito.

Figura 7.12: Mapa de cor das tensões equivalentes de von Mises nos corpos sob carregamentovertical com atrito.

80

7.3 Exemplo 3: Contato de duas sapatas elásticas sob carregamento vertical ehorizontal

Neste exemplo, será considerada a aplicação de força vertical e horizontal. Primeiro aplica-se a carga vertical com as mesmas propriedades do primeiro exemplo e em seguida uma foçatangencial foi aplicada, considerando atrito entre os corpos. A força tangencial considerada é 15

N/mm e o coeficiente de atrito µ = 0.2.

Figura 7.13: Modelagem do problema no programa implementado.

Vale reforçar que os carregamentos vertical e horizontal não podem ser aplicados de umavez. Primeiro a força vertical deve ser aplicada e depois, na outra iteração, a força horizontal. AFigura 7.13 ilustra a modelagem do problema no programa implementado. Na aresta de aplicaçãode carregamento, os nós são fixados na direção y para evitar a rotação da sapata.

81

3 3.5 4 4.5 5 5.5 6 6.5 7

3.5

4

4.5

5

5.5

6

6.5

7

7.5

Figura 7.14: Condições de contorno na aresta superior.

7.3.1 Resultados numéricos

Para esse problema são apresentados resultados para um caso sem atrito e um com atrito noprimeiro passo (carga vertical). A razão para considerar o primeiro passo sem atrito é porque asolução de Mindlin-Cattaneo não considera a tensão de cisalhamento do primeiro passo.

No primeiro caso (considerando sem atrito no passo 1) a tensão cisalhante é a mesma apre-sentada para o primeiro exemplo pois foi considerado o contato sem atrito, e portanto esta tensão énula. Já no segundo passo, onde a carga tangencial é aplicada, o atrito entre os corpos foi conside-rado.

A Figura 7.15 mostra as tensões normal e de cisalhamento no corpo 1 com a malha menosrefinada e compara os resultados numéricos com a solução analítica. A solução analítica para atensão de cisalhamento é a solução de Mindlin-Cattaneo, apresentada anteriormente no capítulo 3.

Na Figura 7.16 se encontram as tensões do corpo 1 com a malha mais refinada. O resultadoobtido mostra ótima concordância entre os resultados numéricos e a solução analítica.

82

Figura 7.15: Tensões normal e cisalhante no corpo 1 com a malha menos refinada e comparaçãocom a solução analítica.

Figura 7.16: Tensões normal e cisalhante no corpo 1 com a malha mais refinada e comparação coma solução analítica.

83

Os deslocamentos na direção x e y para cada corpo estão ilustrados na Figura 7.17. Os deslo-camentos na direção x mostram movimentos dos corpos na direção de aplicação de carregamentotangencial (lado positivo do eixo x, movimento para a direita). A parte de ux em que a curva dosdois corpos coincide mostra a adesão entre eles. Considerando esse movimento como referência, acurva de deslocamento ao longo do eixo x para o corpo 1 é positiva e para o corpo 2 é negativa. Issoé devido a aplicação de força tangencial no corpo 1 que favorece o movimento na direção positivae para corpo 2 é contrária por causa do atrito entre os corpos e a fixação dele, o que faz com que ocorpo resista a este movimento e tente voltar à posição inicial.

Figura 7.17: Deslocamento dos corpos ao longo do eixo x e y no contato sob carregamento verticale horizontal

As curvas de uy também mostram a região de adesão onde as curvas coincidem em tornoda origem do sistema de referência. Os dois corpos têm deslocamento positivo no lado esquerdo enegativo no lado direito. A causa deste deslocamento é o carregamento tangencial no corpo 1, quefaz os corpos se comprimirem do lado direito e se estenderem do lado esquerdo. Um aumento nográfico de uy nos mostra um deslocamento maior no corpo 1, o que é esperado por não ser a regiãode adesão, sendo que este foi o corpo no qual o carregamento foi aplicado.

As tensões equivalentes de von Mises nos corpos na forma de mapa de cor são apresentadosna Figura 7.20.

84

Figura 7.18: Deslocamento dos corpos na direção y no lado esquerdo do corpo.

Figura 7.19: Deslocamento dos corpos na direção y no lado direito do corpo.

85

Figura 7.20: Mapa de cor das tensões equivalente de von Mises nos corpos sob carregamento hori-zontal com atrito.

86

No segundo caso, considerando atrito no passo 1, os resultados foram diferentes. Nesse caso,a força vertical é aplicada gradualmente em 50 incrementos para que o resultado seja mais confiável.Quando a carga externa vertical foi aplicada totalmente, a força horizontal de 2N/mm é aplicada. AFigura 7.21 mostra as tensões e a comparação com a solução analítica. A queda na tensão cisalhanteocorre devido à força tangencial desloca a região de adesão para direita.

−1.5 −1 −0.5 0 0.5 1 1.50

100

200

300

400

500

600

700

x (mm)

Ten

sões

(M

Pa)

σy

τxy

/µ

Figura 7.21: Tensões normal e cisalhante no corpo 1 sob carregamento vertical e horizontal e solu-ção analítica.

A Figura 7.22 mostra os deslocamentos dos corpos. Estes na direção x claramente mostram atransferência da região de adesão para direita, onde as curvas de ux se coincidem. Esta região estárelativamente deslocada para direita em relação ao eixo de referência no meio do corpo. Além disso,o lado direito dos corpos tem um movimento mais próximo ao eixo x enquanto o lado esquerdodo corpo 2 se desloca proporcionalmente menos que o lado esquerdo do corpo 1. Os corpos têmcomportamentos diferentes pois a carga é aplicada diretamente no corpo 1 e o corpo 2 é restrito.Portanto, o corpo 1 se movimenta na direção da aplicação de carga tangencial e o corpo 2 tentaresistir a esse movimento.

As curvas uy apresentam um comportamento parecido com casos de carregamento vertical,ainda que pequenas diferenças possam ser verificadas. Como está apresentado na Figura 7.22, ocomportamento nas extremidades de cada corpo não é simétrico, ou seja, o lado direito dos corpos

87

−8 −6 −4 −2 0 2 4 6 8−0.07

−0.06

−0.05

−0.04

−0.03

−0.02

−0.01

0

0.01

0.02

x (mm)

y (m

m)

uy corpo 1

ux corpo 1

uy corpo 2

ux corpo 2

Figura 7.22: Deslocamento dos corpos ao longo do eixo x e y sob carregamento vertical e horizontalcom atrito.

tem deslocamento maior do que o lado esquerdo. Isso também é devido a carga horizontal ser apli-cada no corpo 1 e a região de adesão entre os corpos. A Figura 7.23 mostra as tensões equivalentede von Mises nos corpos.

88

Figura 7.23: Mapa de cor das tensões equivalente de von Mises nos corpos sob carregamento ver-tical e horizontal com atrito.

89

8 Considerações Finais

8.1 Conclusões

O estudo de contato de dois corpos elásticos em regime de elasticidade infinitesimal (peque-nas deformações) e a implementação de um algoritmo para análise de problemas de contato comuso de métodos dos elementos de contorno foi apresentado. Devido às restrições do problema e apresença de atrito, o problema se torna não linear.

Para as restrições de atrito, a lei de atrito de Coulomb foi considerada. Na presença de atrito,dependendo da força aplicada, a região de contato pode ter zonas de adesão e escorregamento,simultaneamente. Para determinar o estado de contato, um algoritmo baseado no limite de atritofoi proposto. A definição cuidadosa de cada estado tem grande influência na obtenção de tensões edeformações.

Para diminuir a não linearidade do problema, a carga total aplicada é inserida de forma in-cremental. A variação de carga pode mudar o estado de contato definido. Portanto, com cada incre-mento de carga, as restrições são verificadas e os estados são atualizados.

As matrizes são montadas, formando o sistema de equações. Como foi mencionado ante-riormente, esse sistema é um sistema não linear e deve ser resolvido com auxílio de métodos deresolução de sistemas não lineares. Entre os métodos numéricos de minimização conhecidos, o mé-todo de Newton foi implementado e apresentou uma resposta adequada, com número de iteraçõesrelativamente baixo. Este método permite determinar as restrições de contato em cada iteração, demaneira precisa.

A modelagem do problema foi feita com elementos quadráticos contínuos do método dos ele-mentos de contorno (MEC). Para problemas onde o campo de interesse é o contorno, o MEC possuimais vantagens do que qualquer outro método, pois discretiza apenas o contorno do componente.Mesmo que a literatura não recomende o uso de elementos quadráticos por apresentar oscilaçõesna resposta, este tipo de elemento foi utilizado juntamente com o algoritmo proposto. Ao contráriodo mostrado em vários trabalhos da literatura, as respostas obtidas com o uso de elementos de altaordem, no caso desse trabalho, elementos quadráticos, não apresentam oscilações nos resultados

91

das tensões de cisalhamento.

A aplicação de força tangencial foi outro desafio do projeto, pois aumenta mais a complexi-dade do problema. A carga tangencial é aplicada de forma incremental, somente quando a cargavertical atingir seu valor máximo. O método de Newton implementado apresentou pequenas fa-lhas na presença de carga tangencial e precisa de estudos mais detalhados para entender melhor osproblemas relacionados ao método.

O algoritmo apresentado e o programa implementado apresentaram resultados com boa con-cordância com as soluções analíticas disponíveis na literatura. Os resulatados mostram que o MECpode resolver de maneira adequada problemas de contato entre dois corpos. Embora o desenvol-vimento da formulação tenha sido feito para apenas dois corpos, sua extensão para problemasmulti-corpos é direta, com pequenas modificações no código implementado.

8.2 Sugestões para trabalhos futuros

O caso de problemas de contato mecânico é um caso bastante amplo com diversos problemasde interesse na engenharia. Entretanto, o MEC foi pouco explorado nessa área e, portanto, existemmuitos pontos a serem investidos e estudados. Assim, pode-se sugerir os seguintes temas paratrabalhos futuros:

Estudar e estender o programa para caso de aplicação de força tangencial oscilatória e fadigapor fretting;

Extensão do programa para problemas de contato de múltiplos corpos como, por exemplo,cabos condutores;

Extensão para casos tridimensionais;

Extensão para situação de contato dinâmico;

Implementação do MEC isogeométrico e da aproximação cruzada adaptativo (ACA) para aanálise de problemas de contato mecânico.

92

Referências

ABASCAL, R. 2d transient dynamic friction contact problems. i. numerical analysis. Engineeringanalysis with boundary elements, v. 16, n. 3, 227–233, 1995.

ALART, P. e CURNIER, A. A mixed formulation for frictional contact problems prone to newtonlike solution methods. Computer methods in applied mechanics and engineering, v. 92, n. 3,353–375, 1991.

ANDERSSON, T. The boundary element method applied to two-dimensional contact problemswith friction. pp. 239–258, 1981.

ARAUJO, J. e NOWELL, D. The effect of rapidly varying contact stress fields on fretting fatigue.International Journal of Fatigue, v. 24, n. 7, 763–775, 2002.

ARAUJO, J.A. On the initiation and arrest fretting fatigue cracks. 2000. Tese (Doutorado).University of Oxford.

BAIETTO, M.C.; PIERRES, E.; GRAVOUIL, A.; BERTHEL, B.; FOUVRY, S. e TROLLE, B.Fretting fatigue crack growth simulation based on a combined experimental and xfem strategy.International Journal of Fatigue, v. 47, 31–43, 2013.

BARBOSA, R.E. e GHABOUSSI, J. Discrete finite element method for multiple deformable bo-dies. Finite Elements in Analysis and Design, v. 7, n. 2, 145–158, 1990.

BATRA, R. Quasistatic indentation of a rubber-covered roll by a rigid roll. International Journalfor Numerical Methods in Engineering, v. 17, n. 12, 1823–1833, 1981.

93

BRAMHALL, R. Studies in fretting fatigue. 1973. Tese (Doutorado). PhD. Thesis. Universityof Oxford.

BUCHANAN, A.M. e FITZGIBBON, A.W. Damped newton algorithms for matrix factorizationwith missing data. In Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEEComputer Society Conference on, v. 2, pp. 316–322. IEEE, 2005.

CATTANEO, C. Sul contatto di due corpi elastici: distribuzione locale degli sforzi. Rend. Accad.Naz. Lincei, v. 27, n. 6, 342–348, 1938.

CHRISTENSEN, P.W. Algortihms for Frictional Contact Problems Based on MathematicalProgramming. 1997. Tese (Doutorado). Linkoping University, Sweden.

COCU, M. Existence of solutions of signorini problems with friction. International journal ofengineering science, v. 22, n. 5, 567–575, 1984.

CURNIER, A. A theory of friction. International Journal of Solids and Structures, v. 20, n. 7,637–647, 1984.

DUVAUT, G. e LIONS, J.L. Les inéquations en mécanique et en physique. 1972.

FESSLER, H. e OLLERTON, E. Contact stresses in toroids under radial loads. British journal ofapplied physics, v. 8, n. 10, 387, 1957.

FOUVRY, S.; ELLEUCH, K. e SIMEON, G. Prediction of crack nucleation under partial slipfretting conditions. The Journal of Strain Analysis for Engineering Design, v. 37, n. 6, 549–564, 2002.

FOUVRY, S.; KAPSA, P.; SIDOROFF, F. e VINCENT, L. Identification of the characteristic lengthscale for fatigue cracking in fretting contacts. Le Journal de Physique IV, v. 8, n. PR8, Pr8–159,1998.

94

FRANGI, A. e NOVATI, G. Bem–fem coupling for 3d fracture mechanics applications. Compu-tational Mechanics, v. 32, n. 4-6, 415–422, 2003.

FREDHOLM, I. Sur une classe d’équations fonctionnelles. Acta mathematica, v. 27, n. 1, 365–390, 1903.

GIANNAKOPOULOS, A.; LINDLEY, T.; SURESH, S. e CHENUT, C. Similarities of stress con-centrations in contact at round punches and fatigue at notches: implications to fretting fatigue crackinitiation. Fatigue and Fracture of Engineering Materials and Structures, v. 23, n. 7, 561–572,2000.

GINER, E.; SABSABI, M.; RÓDENAS, J.J. e FUENMAYOR, F.J. Direction of crack propagationin a complete contact fretting-fatigue problem. International Journal of Fatigue, v. 58, 172–180,2014.

GONZÁLEZ, J.A.; PARK, K.; FELIPPA, C.A. e ABASCAL, R. A formulation based on localizedlagrange multipliers for bem–fem coupling in contact problems. Computer Methods in AppliedMechanics and Engineering, v. 197, n. 6, 623–640, 2008.

GOODMAN, L. Contact stress analysis of normally loaded rough spheres. Journal of appliedmechanics, v. 29, n. 3, 515–522, 1962.

HARISH, G. e FARRIS, T. Shell modeling of fretting in riveted lap joints. AIAA journal, v. 36,n. 6, 1087–1093, 1998.

HERTZ, H. über die berührung fester elastischer körper (on the contatc of elsatic solids). J. reineund angewandte Mathematik, v. ., ., 1882.

HIBBELER, R. Engineering Mechanics: Statics. Prentice Hall, 2004.

HILLS, D.A. e NOWELL, D. Mechanics of fretting fatique. Springer Netherlands, 1994.

95

HOJJATI TALEMI, R. Numerical modelling techniques for fretting fatigue crack initiationand propagation. 2014. Tese (Doutorado). Ghent University.

IREMAN, P.; KLARBRING, A. e STRÖMBERG, N. Gradient theory of damage coupled to fricti-onal contact and wear, and its numerical treatment. Comput Model Eng Sci, v. 52, n. 2, 125–158,2009.

JOHANSSON, L. Numerical simulation of contact pressure evolution in fretting. Journal oftribology, v. 116, n. 2, 247–254, 1994.

JOHNSON, K.L. Contact mechanics. Cambridge university press, 1987.

KIKUCHI, N. e ODEN, J.T. Contact problems in elasticity: a study of variational inequalitiesand finite element methods, v. 8. siam, 1988.

KLARBRING, A. A mathematical programming approach to three-dimensional contact problemswith friction. Computer Methods in Applied Mechanics and Engineering, v. 58, n. 2, 175–200,1986.

LANDERS, J.A. e TAYLOR, R.L. An augmented alagrangian formulation for the finite elementsolution of contact problems. Relatório técnico, DTIC Document, 1986.

MADGE, J.; LEEN, S.; MCCOLL, I. e SHIPWAY, P. Contact-evolution based prediction of frettingfatigue life: effect of slip amplitude. Wear, v. 262, n. 9, 1159–1170, 2007.

MAN, K.; ALIABADI, M. e ROOKE, D. Bem frictional contact analysis: load incremental tech-nique. Computers & structures, v. 47, n. 6, 893–905, 1993.

MAN, K.W. Contact mechanics using boundary elements. Computational Mechanics, 1994.

MCDOWELL, J.; UHLIG, H.; TIERNEY, W.; MCCLELLAN, A. e HORGER, O.J. Fretting cor-

96

rosion. Anti-Corrosion Methods and Materials, v. 1, n. 7, 247–252, 1954.

MCVEIGH, P.; HARISH, G.; FARRIS, T. e SZOLWINSKI, M. Modeling interfacial conditions innominally flat contacts for application to fretting fatigue of turbine engine components. Internati-onal Journal of Fatigue, v. 21, S157–S165, 1999.

METALS, A. Properties and Selection of Metals. Handbook, 1961.

MINDLIN, R.D. Compliance of elastic bodies in contact. Journal of applied mechanics, v. 16,159–268, 1949.

MOSTOFI, A. e GOHAR, R. Pressure distribution between closely contacting surfaces. Journalof Mechanical Engineering Science, v. 22, n. 5, 251–259, 1980.

MUSKHELISHVILI, N.I. Some basic problems of the mathematical theory of elasticity. 1.Springer Science & Business Media, 1977.

NISHIOKA, K. e HIRAKAWA, K. Fundamental investigations of fretting fatigue:(part 2, frettingfatigue testing machine and some test results). Bulletin of JSME, v. 12, n. 50, 180–187, 1969.

NISHIOKA, K. e KENJI, H. Fundamental investigations of fretting fatigue: Part 6, effects ofcontact pressure and hardness of materials. Bulletin of JSME, v. 15, n. 80, 135–144, 1972.

NISHIOKA, K.; NISHIMURA, S. e HIRAKAWA, K. Fundamental investigations of fretting fati-gue: Part 1, on the relative slip amplitude of press-fitted axle assemblies. Bulletin of JSME, v. 11,n. 45, 437–445, 1968.

NOWELL, D. An analysis of fretting fatigue. 1988. Tese (Doutorado). University of Oxford.

OMBERG, N. An augmented lagrangian method for fretting problems. European Journal ofMechanics-A/Solids, v. 16, 573–593, 1997.

97

PANAGIOTOPOULOS, I.P. A nonlinear programming approach to the unilateral contact-, andfriction-boundary value problem in the theory of elasticity. Ingenieur-Archiv, v. 44, n. 6, 421–432, 1975.

PARIS, F. e GARRIDO, J. An incremental procedure for friction contact problems with the boun-dary element method. Engineering analysis with boundary elements, v. 6, n. 4, 202–213, 1989.

PÁCZELT, I.; KUCHARSKI, S. e MRÓZ, Z. The experimental and numerical analysis of quasi-steady wear processes for a sliding spherical indenter. Wear, v. 274, 127–148, 2012.

PERSSON, A. On the Stress Distribution of Cylindrical Bodies in Elastic Contact. 1964. Tese(Doutorado). PhD Dissertation, Chalmers Teknisha Hogskola.

POPOV, V. Contact mechanics and friction: physical principles and applications. SpringerScience Business Media, 2010.

PUSO, M.A. e LAURSEN, T.A. A mortar segment-to-segment contact method for large deforma-tion solid mechanics. Computer methods in applied mechanics and engineering, v. 193, n. 6,601–629, 2004.

REBEL, G.; PARK, K. e FELIPPA, C. A contact formulation based on localized lagrange mul-tipliers: formulation and application to two-dimensional problems. International Journal forNumerical Methods in Engineering, v. 54, n. 2, 263–297, 2002.

RODRÍGUEZ-TEMBLEQUE, L.; ABASCAL, R. e ALIABADI, M. Anisotropic fretting wearsimulation using the boundary element method. Computer Modeling in Engineering & Scien-ces(CMES), v. 87, n. 2, 127–156, 2012.

RODRÍGUEZ-TEMBLEQUE SOLANO, L. Formulación numérica de la interacción mecánicaentre superficies de sólidos 3d. 2009. Tese (Doutorado). Universidad de Sevilla.

RUIZ, C.; BODDINGTON, P. e CHEN, K. An investigation of fatigue and fretting in a dovetail

98

joint. Experimental Mechanics, v. 24, n. 3, 208–217, 1984.

RUIZ, C. e CHEN, K. Life assessment of dovetail joints between blades and discs in aero-engines.Mechanical Engineering Publications,, v. 15, 187–194, 1986.

RUIZ, C. e NOWELL, D. Designing against fretting fatigue in aeroengines. European StructuralIntegrity Society, v. 26, 73–95, 2000.

SHIGLEY, J.E. Mechanical engineering design. McGraw-Hill, 1972.

SIMO, J.C. e LAURSEN, T. An augmented lagrangian treatment of contact problems involvingfriction. Computers & Structures, v. 42, n. 1, 97–116, 1992.

SINGH, K.P. e PAUL, B. Numerical solution of non-hertzian elastic contact problems. Journal ofApplied mechanics, v. 41, n. 2, 484–490, 1974.

SMIRNOV, V.I. A course of higher mathematics. Pergamon, 1964.

STRÓMBERG, N. A newton method for three-dimensional fretting problems. InternationalJournal of Solids and Structures, v. 36, n. 14, 2075–2090, 1999.

STRÓMBERG, N.; JOHANSSON, L. e KLARBRING, A. Derivation and analysis of a generalizedstandard model for contact, friction and wear. International Journal of Solids and Structures,v. 33, n. 13, 1817–1836, 1996.

TIMOSHENKO, S.P. e GOODIER, J.N. Theory of Elasticity. McGraw-Hill, New York, 1951.

VON MISES, R. Mechanik der festen korper im plastisch-deformablen zustand. Nachrichtenvon der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse,v. 1913, 582–592, 1913.

99

WATERHOUSE, R.B. Fretting fatigue. Elsevier Science & Technology, 1981.

WEI, D.S.; WANG, Y.R. e YANG, X.G. Analysis of failure behaviors of dovetail assemblies dueto high gradient stress under contact loading. Engineering Failure Analysis, v. 18, n. 1, 314–324,2011.

WRIGGERS, P. e LAURSEN, T.A. Computational contact mechanics, v. 30167. Springer, 2006.

WRIGGERS, P.; SIMO, J. e TAYLOR, R. Penalty and augmented lagrangian formulations forcontact problems. Proceedings of the NUMETA, v. 85, 97–105, 1985.

ZAVARISE, G.; WRIGGERS, P. e SCHREFLER, B. On augmented lagrangian algorithms forthermomechanical contact problems with friction. International Journal for Numerical Methodsin Engineering, v. 38, n. 17, 2929–2949, 1995.

100

A Código computacional da solução analítica Mindlin-Cattaneo

% Soluções analíticas para o problema de contato entre um cilindro e uma

% base rígida

% Solução para a tensão normal (p): Solução de Hertz

% Solução para a tensão de cisalhamento (q): Solução de Mindlin-Cattaneo

clear

close all

cargav=100; % Valor da carga aplicada na direção vertical

cargah=15; % Valor da carga aplicada na direção horizontal

w=6.5; % Espessura da sapata

R=70; % Raio da sapata

P=2*w*cargav; % Carga total aplicada na direção vertical

Q=2*w*cargah; % Carga total aplicada na direção horizontal

ni=0.33; % Coeficiente de Poisson

E=73.4e3; % Módulo de elasticidade

East=E/(1-ni); % E no estado plano de deformação

a=sqrt(4*R*P/East/pi); % Comprimento da região de contato

xno=-a:a/100:a; % Coordenadas x dos nós

po=2*P/(pi*a); % Valor máximo da tensão normal na reginão de contato

qo=Q/(2*pi*a^2);

mi=cargah/cargav; % Valor mínimo possível para o coeficiente de atrito

for j=1:4 % Loop sob os coeficientes de atrito

c=a*sqrt(1-abs(Q/(mi*P))); % Comprimento da região em adesão

vetmi(j)=mi; % Vetor que guarda os valores dos coeficientes de atrito

for i = 1:size(xno,2) % Loop sobre os nós da região de contato

if(abs(xno(i))<=a) % Verifica se o nó está em contato

p(i)=po*sqrt(1-(xno(i)/a)^2); % tensão normal

q(i,j)=mi*p(i); % tensão de cisalhamento (região em slip)

if(abs(xno(i))<=c) % Verifica se o nó está na região de adesão

qlinha=mi*po*c/a*sqrt(1-(xno(i)/c)^2);

q(i,j)=q(i,j)-qlinha; % tensão de cisalhamento (região

% em stick)

end

else

p(i)=0; % Valor da tensão normal na região que não está em

% contato

q(i,j)=0; % Valor da de cisalhamento normal na região que

101

% não está em contato

end

end

mi=mi+0.05; % Incremento do coeficiente de atrito

end

figure

plot(xno/a,p/po,'bd-',xno/a,q(:,1)/(vetmi(1)*po),'rs-', ...

xno/a,q(:,2)/(vetmi(2)*po),'ko-',xno/a,q(:,3)/(vetmi(3)*po),'g*-', ...

xno/a,q(:,4)/(vetmi(4)*po),'mp-')

leg1=['q para Q/(\mu P) =',num2str(Q/(vetmi(1)*P))];

leg2=['q para Q/(\mu P) =',num2str(Q/(vetmi(2)*P))];

leg3=['q para Q/(\mu P) =',num2str(Q/(vetmi(3)*P))];

leg4=['q para Q/(\mu P) =',num2str(Q/(vetmi(4)*P))];

legend('p/p_o',leg1,leg2,leg3,leg4);

ylabel('Tensões normalizadas')

xlabel('\it x/a')

title('Tensões normal e de cisalhamento para diferentes situações')

102