Escola Politécnica da Universidade de São Paulo Departamento de Energia e Automação Elétrica

O MÉTODO DOS ELEMENTOS FINITOS COMO FERRAMENTA PARA

PROJETO DE MALHAS DE ATERRAMENTO DE SUBESTAÇÕES

DANILO AUGUSTO SALGADO MARCELO PEREIRA PINTO

São Paulo Junho de 2009

Danilo Augusto Salgado Marcelo Pereira Pinto

O MÉTODO DOS ELEMENTOS FINITOS COMO FERRAMENTA PARA

PROJETO DE MALHAS DE ATERRAMENTO DE SUBESTAÇÕES

Trabalho de formatura apresentado à Escola Politécnica da Universidade de São Paulo para a obtenção do grau de graduado em Engenharia Elétrica. Área de Concentração: Sistemas de Potência Orientadores: Prof. Dr. José Roberto Cardoso Profª. Drª. Viviane Cristine Silva

São Paulo 2009

Danilo Augusto Salgado Marcelo Pereira Pinto

O MÉTODO DOS ELEMENTOS FINITOS COMO FERRAMENTA PARA

PROJETO DE MALHAS DE ATERRAMENTO DE SUBESTAÇÕES

Trabalho de formatura apresentado à Escola Politécnica da Universidade de São Paulo para a obtenção do grau de graduado em Engenharia Elétrica. Área de Concentração: Sistemas de Potência Orientadores: Prof. Dr. José Roberto Cardoso Profª. Drª. Viviane Cristine Silva

São Paulo 2009

DANILO AUGUSTO SALGADO MARCELO PEREIRA PINTO

O MÉTODO DOS ELEMENTOS FINITOS COMO FERRAMENTA PARA

PROJETO DE MALHAS DE ATERRAMENTO DE SUBESTAÇÕES

_________________________________________________ PROF. DR. JOSÉ ROBERTO CARDOSO

_________________________________________________ PROF. DRA. VIVIANE CRISTINE SILVA

FICHA CATOLOGRÁFICA SALGADO, DANILO AUGUSTO PINTO, MARCELO PEREIRA

CÁLCULO DE SISTEMAS DE ATERRAMENTO ELÉTRICO ATRAVÉS DE ELEMENTOS FINITOS / D. A. SALGADO E M. P. PINTO – SÃO PAULO, 2009.

TRABALHO DE FORMATURA (GRADUAÇÃO) – ESCOLA POLITÉCNICA DA UNIVERSIDADE DE

SÃO PAULO. DEPARTAMENTO DE ENGENHARIA DE ENERGIA E AUTOMAÇÃO ELÉTRICAS. 1. ELETROMAGNETISMO 2. MÉTODO DOS ELEMENTOS FINITOS 3. ATERRAMENTO 4.

SISTEMAS ELÉTRICOS DE POTÊNCIA I. UNIVERSIDADE DE SÃO PAULO. ESCOLA

POLITÉCNICA. DEPARTAMENTO DE ENERGIA E AUTOMAÇÃO ELÉTRICAS.

i

AGRADECIMENTOS

Ao Prof. Dr. José Roberto Cardoso, orientador deste trabalho, professor de didática

inigualável, pela amizade e por ter demonstrado a centenas de alunos, inclusive os autores,

a grandeza do eletromagnetismo de maneira clara e por nos lembrar constantemente da

nobreza da profissão do engenheiro.

À Profª. Drª. Viviane Cristine Silva, orientadora deste trabalho, pela atenção e dedicação,

por fornecer os subsídios para a construção desta pesquisa e pela paciência.

Aos professores da graduação do Departamento de Energia e Automação Elétricas, pelo

esforço em manter a excelência do ensino na Escola Politécnica. Estas pessoas foram

fundamentais na formação de centenas de engenheiros que hoje contribuem para o

desenvolvimento do nosso país. Apesar de todas as dificuldades, seus esforços não foram

em vão.

Aos ex-alunos e professores que contribuíram pelas melhorias realizadas neste

departamento nos últimos anos. O trabalho realizado por eles tornou este departamento

referência dentro de nossa escola.

Às amizades que construímos dentro da escola, muitas das quais serão carregadas por toda

a vida. Estiveram sempre presentes nos momentos felizes e nos apoiaram nos momentos de

dificuldade.

Aos amigos da Keep Flying, uma equipe de competição da Escola Politécnica que os

autores deste trabalho participaram ativamente. Na Keep Flying fizemos grandes amigos e

aprendemos a aplicar a engenharia na prática, tendo a oportunidade de levar o nome da

nossa escola ao topo ao alçar o primeiro lugar geral em uma competição de engenharia

aeronáutica de âmbito mundial promovida nos Estados Unidos da América.

Aos nossos avôs, que não tiveram a oportunidade de estudar, mas que deram condições a

seus filhos e netos para que o tivessem, apesar da vida árdua. A eles somos eternamente

gratos.

ii

DEDICATÓRIA

Às nossas famílias e

às nossas namoradas

iii

“If I have seen further it is by standing on the shoulders of giants”

Isaac Newton

iv

RESUMO

Este trabalho tem como objetivo demonstrar aplicações práticas do método numérico de

elementos finitos no projeto e análise de sistemas reais de aterramento, em regime

permanente. A principal preocupação é comparar os resultados numéricos obtidos com os

cálculos analíticos e semi-empíricos e com as medições feitas em campo para verificar a

validade e aplicabilidade do método dos elementos finitos para sistemas de aterramento.

O trabalho não se resume somente ao método dos elementos finitos, mas também

discorre brevemente sobre o solo e sua resistividade, os sistemas de aterramento mais

comuns e seus cálculos analíticos ou semi-empíricos através da norma mais utilizada no

setor (IEEE Std 80). Essa introdução ao tema do aterramento se fez necessária para que as

posteriores comparações pudessem ser feitas.

O método dos elementos finitos é descrito em detalhes, desde sua formulação até os

aspectos práticos de implantação, como geração de geometrias e malhas para simulação.

Para essas etapas citadas, utilizou-se um aplicativo livremente distribuído, o GMSH,

ferramenta que facilitou imenso o desenrolar deste trabalho.

Utilizou-se também uma técnica moderna para minimizar os erros devido ao

truncamento do domínio, o PML, e para as simulações utilizou-se o solver do GROUND

3D do laboratório de eletromagnetismo aplicado (LMAG).

Os autores se utilizaram de dados de subestações elétricas reais, tanto dados de projeto,

como também experimentais, obtidos em medições em campo para realizar as comparações

entre os métodos de projeto. Ao final, discute-se a validade e a aplicabilidade do método

dos elementos finitos no projeto de aterramento de subestações elétricas.

Palavras chaves: Eletromagnetismo. Método dos Elementos Finitos. Aterramento.

Sistemas Elétricos de Potência.

v

ABSTRACT

The goal of this work is to demonstrate practical applications of the finite element method

applied to the design and analysis of real grounding systems in steady-state conditions. The

numerical results were compared to the analytical and semi-empirical methods and the in-

locus measurements in order to study the feasibility and applicability of the finite element

method to grounding systems.

This work also addresses the soil and its resistivity, the most simple grounding techniques

and its formulation, analytical, semi-empirical and accordingly to the IEEE Std 80

Standard. This brief introduction is necessary due to the comparisons done in the further

chapters.

The finite element method is explained in details, from its formulation to its practical

applications, for example, the geometry and mesh generation. These steps were

accomplished with free software called GMSH.

A modern technique (PML, Perfectly Matched Layer) was employed to minimize the errors

due to the domain truncation. The simulations were carried with the GROUND 3D solver

from the LMAG Laboratory.

The authors used data from real substations, from design data to in-locus measurements, in

order to compare the design methods described in this work. At the end, the feasibility of

the finite element method applied to grounding systems is discussed.

Keywords: Electromagnetism, Finite Element Method, Grounding, Power Systems.

vi

SUMÁRIO

1. Introdução................................................................................................................................................... 1 2. Solo............................................................................................................................................................. 5

2.1 A resistividade do solo......................................................................................................................... 5 2.2 Medida de Resistividade e Estratificação do Solo ............................................................................... 9

3. Métodos Analíticos................................................................................................................................... 12 3.1 Resistência de Hastes Verticais.......................................................................................................... 13 3.2 Resistência de condutores dispostos horizontalmente........................................................................ 15 3.3 Resistência de uma malha de aterramento segundo a IEEE Std 80.................................................... 18

4. Método dos Elementos Finitos ................................................................................................................. 20 4.1 Métodos Numéricos para Eletromagnetismo ..................................................................................... 20 4.2 Aplicação do MEF em Sistemas de Aterramento............................................................................... 24

4.2.1 Formulação Matemática ........................................................................................................... 24 4.2.2 Divisão do Domínio de Estudo................................................................................................. 27 4.2.3 Funções Interpoladoras e Funções de Forma............................................................................ 29 4.2.4 Formulação Matricial do Problema .......................................................................................... 31 4.2.5 A Matriz Global e a Potência Dissipada................................................................................... 33 4.2.6 A Solução do Sistema............................................................................................................... 34 4.2.7 Resistência do Sistema de Aterramento.................................................................................... 35

4.3 PML (Perfectly Matched Layers)....................................................................................................... 36 4.4 GMSH................................................................................................................................................ 40

5. Aplicações e Resultados ........................................................................................................................... 44 5.1 Comparação com Métodos Analíticos ............................................................................................... 45 5.2 Comparação com o Método da Norma IEEE Std 80.......................................................................... 49 5.3 Comparação com uma Malha Real .................................................................................................... 52 5.4 Resumo das Comparações.................................................................................................................. 55

6. Conclusões................................................................................................................................................ 57 7. Bibliografia............................................................................................................................................... 61

vii

LISTA DE ILUSTRAÇÕES

Ilustração 2-1 – principais grupos minerais – fonte: Wikipedia (www.wikipedia.org)...................................... 6 Ilustração 2-2 – Método de Wenner ................................................................................................................... 9 Ilustração 3-1 – distribuição das equipotenciais em uma haste simples ........................................................... 13 Ilustração 3-2 – condutor disposto horizontalmente......................................................................................... 15 Ilustração 3-3 condutores dispostos em L e em cruz........................................................................................ 16 Ilustração 3-4 – parâmetros geométricos de uma malha retangular.................................................................. 18 Ilustração 4-1 domínio de estudo truncado....................................................................................................... 25 Ilustração 4-2 elementos típicos para o método dos elementos finitos............................................................. 28 Ilustração 4-3 truncamento do domínio por material absorvedor anisotrópico PML ....................................... 37 Ilustração 4-4 – geração da geometria.............................................................................................................. 41 Ilustração 4-5 – geração da malha .................................................................................................................... 41 Ilustração 4-6 – adensamento de volumes nos pólos e no entreferro ............................................................... 41 Ilustração 4-7 – distribuição de pressão nas palhetas de um estágio de uma turbina para geração termoelétrica..................................................................................................................................................... 42 Ilustração 4-8 – geometria criada no Gmsh e sua malha de elementos finitos ................................................. 43 Ilustração 5-1 – distribuição do potencial elétrico para o caso da haste delgada de 32 m, ρ = 450 Ω.m e r = 4 mm.................................................................................................................................................................... 46 Ilustração 5-2 – distribuição de potencial elétrico para quarto casos que possuem solução analítica .............. 48 Ilustração 5-3 – geometria da malha utilizada na comparação do MEF com a norma IEEE Std 80 ............... 49 Ilustração 5-4– distribuição de potencial para os seis perfis de solo para a malha 16 x 16 m.......................... 51 Ilustração 5-5 - geometria da malha de aterramento simulada neste trabalho .................................................. 52 Ilustração 5-6 – distribuição de potencial elétrico na malha............................................................................. 54 Ilustração 5-7 – distribuição de potencial ao longo de uma diagonal da malha ............................................... 54

LISTA DE TABELAS

Tabela 2-1 – Resistividade para alguns tipos de solo ......................................................................................... 7 Tabela 2-2 – resistividade de um solo arenoso para diversos índices de umidade ............................................. 7 Tabela 2-3 – resitividade de um solo arenoso úmido para diversas temperaturas .............................................. 7 Tabela 3-1 – comparação entre geometrias ...................................................................................................... 17 Tabela 4-1– Equações para o eletromagnetismo geral e eletrocinética nas suas formas diferenciais............... 24 Tabela 5-1 comparação de métodos para uma haste de 32 m de comprimento, ρa = 450 Ω.m e r = 4 mm...... 45 Tabela 5-2 – comparação de métodos para quatro casos analíticos.................................................................. 48 Tabela 5-3 – comparação do MEF com a norma IEEE Std 80 para a malha 16 x 16 m................................... 50 Tabela 5-4 – comparação de métodos para a malha de uma subestação típica, ρa = 450 Ω.m e r = 4 mm....... 53

viii

LISTA DE SÍMBOLOS

a parâmetro do PML Rh resistência de terra da haste [Ω]

Am área ocupada pela malha [m²] Rm resistência de terra da malha [Ω] B densidade de fluxo magnético [Wb/m²] ou [T] r raio do condutor [m]

d distância até a haste [m] s parâmetro de “esticamento” das coordenadas

D deslocamento [C/m²] t tempo [s]

E campo elétrico [V/m] Vei potencial no nó i do elemento e G matriz global V potencial elétrico [V]

Gij coeficiente da matriz global x,y,z coordenadas cartesianas

H campo magnético [A/m] ρa resistividade aparente do solo [Ω.m]

hm profundidade da malha [m] ρ1 resistividade da primeira camada [Ω.m]

I corrente de defeito [A] ρ2 resistividade da segunda camada [Ω.m]

J densidade de corrente de condução [A/m²] ρ densidade volumétrica de carga [C/m³]

L comprimento da haste [m] σ condutividade elétrica [S/m]

Lt comprimento total de cabos da malha [m] Ω domínio de estudo

n número de nós do elemento ГD fronteira sujeita à condição de Dirichlet

N número de elementos da malha ГN fronteira sujeita à condição de von Neumann

n vetor normal à superfície αi função de forma do elemento

p profundidade [m] σ tensor de condutividade elétrica [S/m]

Pe potência dissipada no elemento e [W] Λ matriz de transformação do domínio

pe densidade de potência dissipada no elemento e [W/m³]

Índices: e - elemento

ABREVIATURAS GMSH programa gerador de malhas de elementos finitos

ICCG Incomplete Cholesky Conjugated Gradients – Método dos Gradientes Conjugados por

Decomposição Incompleta de Cholesky

LMAG Laboratório de Eletromagnetismo Aplicado Prof. Dr. Áurio Gilberto Falcone

MEF Método dos Elementos Finitos

MDF Método das Diferenças Finitas

MRT Monofásico com Retorno pela Terra

PML Perfect Matched Layers

TLM Transmission Line Method (Método das Linhas de Transmissão)

1

1. INTRODUÇÃO

Os sistemas elétricos de potência têm por principal função conduzir largas quantidades

de energia das usinas de geração até os grandes centros consumidores sob a forma de

eletricidade. Em razão da resistividade dos condutores que transportam essa energia, há

perdas no sistema de transmissão e distribuição e para reduzi-las os níveis de tensão elétrica

que os sistemas operam são elevadíssimos.

Os níveis de tensão, mesmo os mais baixos como o do sistema de distribuição,

representam uma ameaça à vida humana e essa é a principal razão da existência dos

sistemas de aterramento. Tais sistemas permitem a condução adequada das correntes de

falta ao solo, fazendo com que os seres humanos fiquem sujeitos a distribuições de

potenciais elétricos aceitáveis, que não lhes causem danos. Porém, essa não é a única razão

que justifica a necessidade do aterramento. Há muitos outros motivos que justificam sua

necessidade, como a proteção de equipamentos elétricos, conduzir as descargas

atmosféricas para a terra, prover caminho para as correntes de falta, prover um caminho de

retorno de corrente nos sistemas MRT e, também, melhorar a sensibilidade dos

equipamentos de proteção de forma a melhor proteger os circuitos.

Há tempos, os projetos das malhas de aterramento têm sido calculados de acordo com

normas específicas, utilizando-se de modelos bem simplificados do solo, tabelas e equações

semi-empíricas que foram compiladas ao longo dos anos. Tais modelos não prevêem a real

distribuição de potencial ao longo do solo, não lidam com geometrias complexas e por

muitas vezes, em virtude das simplificações assumidas, provêem resultados de acuracidade

duvidosa. Contudo, isso não significa, necessariamente, que esses métodos não sejam

aplicáveis, porém há de se reconhecer suas severas limitações e buscar outros métodos que

a complementem. Surge então a idéia do método numérico.

O vertiginoso crescimento da capacidade computacional nas últimas décadas

possibilitou a resolução numérica de muitos problemas práticos de engenharia, e não foi

diferente com os sistemas de aterramento. Há cerca de duas décadas, o Laboratório de

2

Eletromagnetismo Aplicado (LMAG) da Escola Politécnica da USP, deu um grande passo

ao criar ferramentas computacionais para a resolução de problemas eletromagnéticos. Em

particular, esse laboratório foi pioneiro na aplicação de métodos numéricos para sistemas

de aterramento, tendo contribuído com inúmeros artigos de relevância internacional.

O método numérico permite que qualquer geometria seja tratada, simples ou complexa,

desde que não se extrapole o limite da memória do computador. Essa limitação ocorre em

virtude da natureza do método, pois é impossível retratar computacionalmente um domínio

de estudo infinito. Adicionalmente, o método numérico resolve os campos elétricos e

magnéticos em inúmeros pontos do domínio de estudo, o que é de especial interesse para os

problemas de aterramento. Nesses sistemas, a distribuição de potencial elétrico na

superfície do solo determina a máxima diferença de potencial que estaria sujeita uma

pessoa na subestação quando da ocorrência de uma falta a terra. Por exemplo, os potenciais

de toque e de passo são parâmetros de segurança importantíssimos que tomam parte no

projeto de uma subestação. Essas duas vantagens tornam atrativo o emprego do método

numérico na análise e projeto de um sistema de aterramento.

Ainda que os métodos numéricos sejam mais generalistas, eles possuem limitações que

restringem a sua aplicação indiscriminada. A capacidade de memória talvez seja a maior

limitação do modelo numérico, pois quanto mais pontos forem analisados dentro de um

domínio, menores serão os erros de discretização. Além disso, a memória limita o tamanho

do domínio de estudo, causando erros de truncamento. Outra grande limitação é a

capacidade de processamento, pois as equações são geralmente muito complexas e tomam

muito tempo para serem resolvidas. Há ainda o fato da maioria das rotinas serem iterativas,

o que aumenta ainda mais o tempo de processamento. Não é de se estranhar que, por muitas

vezes, se leve muitas horas até uma simulação ser concluída.

Alguns métodos numéricos, como o de elementos finitos, exigem que o domínio de

estudo seja dividido em pequenos volumes, dentro dos quais serão resolvidas as equações

diferenciais que regem o modelo. O processo dessa divisão, dita criação da malha de

elementos finitos, toma um tempo considerável, mesmo que gerada automaticamente. As

características dessa malha podem influenciar significativamente os resultados, a ponto de

uma malha ruim gerar resultados inutilizáveis. Este é outro ponto importantíssimo a ser

3

levado em consideração: a ferramenta numérica por si só não é eficiente na solução dos

problemas, cabe ao usuário saber utilizá-la adequadamente e também cabe a ele interpretar

os resultados e julgar se são aplicáveis ou não.

Nestes últimos parágrafos pretendeu-se mostrar que não há método e nem modelo

perfeito, cada qual com suas vantagens e limitações. Apesar deste trabalho se concentrar

nos modelos numéricos, em momento algum se pode afirmar que eles são superiores aos

métodos analíticos e experimentais. A justificativa deste trabalho é demonstrar que o

método numérico é aplicável e conveniente a síntese e a análise de sistemas de aterramento,

principalmente aqueles com geometrias complexas.

Neste trabalho, nos próximos capítulos, também se demonstra que há outros parâmetros

importantíssimos em um bom projeto, como a determinação da faixa de valores que a

resistividade do solo pode assumir. Pode-se até dizer que de nada servem os melhores

métodos sem os melhores modelos. O solo é complexo e seu modelo é por demasiado

simplificado e, apesar de não ser o foco deste trabalho, um capítulo inteiro foi dedicado ao

seu estudo, sua modelagem e a sua medição.

O capítulo 3 trata muito brevemente dos métodos analíticos e semi-empíricos que são

utilizados nos projetos de aterramento, para que o método numérico possa ser validado. De

certa forma, todo modelo numérico deve ser calibrado, testado e comparado com modelos

exatos, para que se tenha o mínimo de confiança nos resultados. As equações analíticas

serão listadas sem alguma demonstração formal, pois não faz parte do escopo.

O capítulo 4 discorre sobre os quatro principais métodos numéricos utilizados na

solução de problemas eletromagnéticos e demonstra-se, brevemente e sem muita

formalidade, que o método dos elementos finitos (MEF) é o mais adequado para o estudo

do aterramento em regime estacionário, representativo para baixas freqüências. Nota-se que

quase todos os estudos numéricos de sistemas de aterramento elétrico se utilizam do MEF,

salvo exceções de estudos em altas-freqüências ou transientes no domínio do tempo.

Depois, retoma-se o método dos elementos finitos, demonstrando sua fundamentação

matemática e sua implantação sob a forma de rotinas de computador. Além disso, comenta-

4

se também da otimização do processamento, tomando partido da esparsidade da matriz

global. Então se aborda a aplicação do PML – Perfect Matched Layer – um artifício

matemático que reduz drasticamente os erros devido ao truncamento do domínio. Discute-

se também a utilização do GMSH, aplicativo livre para geração de malhas, pré e pós-

processamento dos resultados.

No capítulo 5 o MEF é validado com modelos analíticos e então aplicado a malhas de

aterramento reais de subestações de grandes indústrias e usinas hidroelétricas. Comparam-

se então os resultados obtidos numericamente com os cálculos analíticos e semi-empíricos.

Neste capítulo procura-se demonstrar a validade do método numérico na síntese e na

análise de sistemas de aterramento.

Por fim, o capítulo 6 traz as conclusões e comentários finais sobre a aplicação do MEF,

resume os resultados obtidos, indica possíveis complementações a este trabalho, discorre

sobre as dificuldades encontradas na implementação do método e pondera sobre a possível

utilização de softwares de otimização acoplados com o MEF para geração de malhas de

aterramento ótimas.

5

2. SOLO

2.1 A RESISTIVIDADE DO SOLO

A estrutura do solo afeta significativamente um sistema de aterramento elétrico. Além

disso, sua composição varia muito com o tempo e com o clima, fazendo com que os

parâmetros elétricos do aterramento se modifiquem. Um bom projeto deve levar essas

questões em consideração, para que todos os critérios de segurança sejam satisfeitos.

O tipo de solo, a umidade e a temperatura podem fazer com que sua resistividade

assuma valores tão baixos quanto alguns ohms-metro e cresça até valores tão altos quanto

dezenas de milhares de ohms-metro. Essa constante incerteza na resistividade do solo é um

dos grandes problemas no projeto de um sistema de aterramento elétrico e, mesmo através

de ensaios em campo, não se assegura que ela se manterá estável com o decorrer do tempo.

Soma-se a isso a natureza complexa do estudo dos campos eletromagnéticos e a ausência de

modelos bem representativos da condução elétrica no solo e tem-se então uma noção da

dificuldade de se projetar uma malha de aterramento. Portanto, para compreender como o

solo afeta o projeto de um sistema de aterramento elétrico, precisa-se compreender o que

vem a ser o solo em si. Neste capítulo, busca-se entender como é sua estrutura e

composição e, principalmente, entender os mecanismos de condução elétrica através dele.

Para as aplicações de engenharia elétrica, entende-se por solo como tudo o que está

abaixo do plano da terra, vulgarmente, tudo o que está debaixo de nossos pés. Porém, em

geologia, o solo se refere apenas à camada de material que se situa acima das rochas

consolidadas. Essas rochas surgiram a milhões de anos, e a ação da chuva e do vento da

atmosfera primitiva as lixiviaram causando um acúmulo de material, formando assim um

solo primitivo argiloso. O processo natural de lixiviação, decomposição, acumulação e

sedimentação fizeram com que a profundidade do solo aumentasse com o passar do tempo.

Formaram-se diversas camadas, cada qual com suas características químicas e físicas, como

a resistividade elétrica.

6

O solo contém três grupos minerais: a argila, que é composta por grãos muito finos e

que consegue reter muita água, o silte (ou limo) que é composto por grãos menos finos e

que retém pouca água e a areia, que são grãos espessos e que retém quase nada de água.

Como a água é o principal responsável pela condução elétrica no solo, é de se esperar que a

resistividade de um solo argiloso seja significativamente menor do que a de um solo siltoso

ou arenoso. Pedregulhos, seixos, cascalhos e brita, por exemplo, não conseguem reter água

alguma e por isso são utilizados em subestações para isolar as pessoas do contato direto

com o solo mais úmido. Além disso, as rochas consolidadas são praticamente impermeáveis

e conduzem quase nada. A Ilustração 2-1 mostra os grupos mineirais. A

Tabela 2-1 provê resistividades representativas para alguns tipos de solos com suas

umidades típicas.

Argila

Silte

Areia

Brita

Ilustração 2-1 – principais grupos minerais – fonte: Wikipedia (www.wikipedia.org)

7

Tabela 2-1 – Resistividade para alguns tipos de solo

Tipo de Solo Resistividade (Ω.m) Lama 5 a 100

Terra de jardim com 50% de umidade 140 Terra de jardim com 20% de umidade 480

Argila seca 1 500 a 5 000 Argila com 40% de umidade 80 Argila com 20% de umidade 330

Areia seca 3 000 a 5 000 Calcário compacto 1 000 a 5 000

Granito 1 500 a 10 000

Fonte: Kindermann e Campagnolo [1]

A água presente no solo favorece a condução da corrente elétrica, pois os sais oriundos

das rochas e presentes no solo acabam por se dissolver, formando íons, um meio eletrolítico

favorável à condução da corrente iônica. Portanto é de se esperar que quanto maior a

umidade presente no solo, maior será a sua capacidade de condução de corrente elétrica. De

fato, qualquer tipo de solo conduz corrente alguma quando perfeitamente seco, porém

conduz razoavelmente quando um mínimo de umidade está presente.

A temperatura também influi significativamente na resistividade do solo, pois determina

a dispersão das ligações iônicas entre os grânulos de terra no solo. Quando a água congela a

resistividade cresce vertiginosamente com o decréscimo de temperatura. Para se ter uma

noção do efeito da umidade e da temperatura, toma-se como exemplo o solo arenoso, como

os dados apresentados na Tabela 2-2 e Tabela 2-3.

Tabela 2-2 – resistividade de um solo arenoso para diversos índices de umidade

Umidade (% do peso)

Resistividade (Ω.m) de um solo arenoso

0,0 10 000 000 2,5 1 500 5,0 430 10,0 185 15,0 105 20,0 63 30,0 42

Fonte: Kindermann e Campagnolo [1]

Tabela 2-3 – resitividade de um solo arenoso úmido para diversas temperaturas

Temperatura (ºC)

Resistividade (Ω.m) de um solo arenoso

20 72 10 99

0 (água) 138 0 (gelo) 300

-5 790 -15 3 300

Fonte: Kindermann e Campagnolo [1]

8

As tabelas Tabela 2-2 e Tabela 2-3 podem sugerir que a

Tabela 2-1 não tem utilidade alguma na previsão da resistividade do solo a ser estudado.

Em verdade, cada tipo de solo tem uma faixa de umidade e consequentemente uma faixa de

valores possíveis para a sua resitividade. Por exemplo, a argila é tão fina que naturalmente

sempre haverá um mínimo de umidade, ainda mais por estar localizada mais

profundamente, onde a variação de umidade e temperatura é menor, e sempre haverá um

máximo de umidade, pois além de certa quantidade de água esse solo saturaria. E isso

ocorre com todos os horizontes mais profundos. Porém, a concentração de água (umidade)

dos solos mais superficiais e sua temperatura variam significativamente com o clima, as

chuvas e o período do ano. Se o solo superficial tiver baixa capacidade de retenção de água,

como um solo arenoso, essas variações serão ainda maiores.

Esse capítulo procura demonstrar que sempre haverá uma grande incerteza quanto à

resistividade do solo e que o modelo de solo homogêneo, perfeitamente estratificado e

contínuo não tem correspondência real. Por melhor que seja o método aplicado a análise de

um sistema de aterramento, seja ele numérico ou analítico, ele jamais representará

fielmente o fenômeno físico e carregará consigo grandes incertezas que devem ser

enfatizadas. Contudo, até o presente momento não há outras modelos para o solo, pelo

menos que sejam implantáveis, e o modelo de camadas homogêneas, contínuas e

perfeitamente estratificadas satisfazem a necessidade dos projetos de aterramento, dando-

lhos uma visão qualitativa e, por vezes, quantitativa com ressalvas.

9

2.2 MEDIDA DE RESISTIVIDADE E ESTRATIFICAÇÃO DO SOLO

Conhecer a resistividade do solo para várias profundidades é extremamente importante

para um bom projeto de um sistema de aterramento. As medições devem ser feitas na

região a ser construído o sistema, ao longo de várias direções e tomando-se várias medidas

de resistência, para obter uma média representativa e com o menor desvio possível.

Existem três experimentos largamente divulgados para se medir a resistividade do solo para

diversas profundidades: o método de Lee, o de Schlumberger-Palmer e o de Wenner. Este

último é sem dúvida alguma o mais utilizado, por ser muito prático e simples.

Para o método de Wenner, são necessárias quatro hastes com cerca de meio metro de

comprimento e de 10 a 15 mm de diâmetro, de material condutor, isenta de gordura ou

oxidação que possa atrapalhar o bom contato elétrico com o solo. Estas hastes devem ser

dispostas em linha reta ao longo do solo a ser estudado, enterradas à mesma profundidade,

entre 20 e 30 cm, e espaçadas de uma distância a, conforme mostrado na Ilustração 2-2.

Ilustração 2-2 – Método de Wenner

10

Injetando-se uma corrente através das hastes externas, mede-se a diferença de potencial

entre as duas hastes internas, e a partir da relação tensão por corrente injetada surge uma

resistência elétrica que é proporcional à resistividade do solo. Existem aparelhos próprios

para executar esse procedimento, os terrômetros, sendo o mais famoso deles o Megger.

A referência [1], demonstram a fórmula de Palmer para um solo homogêneo (de

resistividade constante), que relaciona os parâmetros geométricos e a resistência medida

com a resistividade do solo. A mesma referência simplifica a fórmula de Palmer para casos

em que o espaçamento é muito grande frente ao comprimento da hastes, que resulta na

Equação 2-1.

Raπρ 2= Equação 2-1

O método de Wenner não contempla apenas o caso de solo homogêneo. Na prática, solo

algum é dessa forma. Segundo [1], a resistividade obtida pelo método de Wenner é

representativa da resistividade do solo a uma profundidade igual ao espaçamento entre

hastes:

O método [de Wenner] considera que praticamente 58% da distribuição de corrente que passa entre as hastes externas ocorre a uma profundidade igual ao espaçamento entre as hastes... A corrente atinge uma profundidade maior, com uma área de dispersão grande, tendo, em conseqüência, um efeito que pode ser desconsiderado. Portanto, para efeito do Método de Wenner, considera-se que o valor da resistência elétrica lida no aparelho [terrômetro] é relativa a uma profundidade a do solo.

(Ref. [1], 1995, p. 21)

Para efeito de modelagem do sistema de aterramento, tanto analítica quanto numérica,

uma distribuição contínua e variável de resistividade em função da profundidade, como a

obtida através do método de Wenner, não se mostra prática. Costuma-se, então, a dividir o

solo em camadas de resistividade constante, representativas de certa faixa de profundidade.

O solo pode ser estratificado em quantas camadas forem necessárias, em geral não mais que

duas ou três, pois as camadas muito mais profundas que o sistema de aterramento tem

pouca influência sobre a resistividade aparente. A referência [1] demonstra vários métodos

de estratificação, os de duas camadas usando curvas, usando técnicas de otimização,

métodos de camada dupla simplificados, o método de Yokogawa e métodos para várias

camadas como o de Pirson. Não é escopo desse trabalho demonstrar algum método de

11

estratificação, basta apenas o leitor estar ciente de que é possível, adequado e recomendável

representar o solo em camadas de resistividade constante.

Em suma, conclui-se que a resistividade do solo depende de sua composição, da

temperatura e principalmente de sua umidade. Além disso, mesmo que a resistividade varie

ao longo da profundidade, é adequado modelar o solo em duas ou três camadas

representativas das profundidades de interesse.

12

3. MÉTODOS ANALÍTICOS

É muito comum na ciência construírem-se modelos complexos a partir de modelos

muito mais simples, bem conhecidos e validados. Com o advento da análise numérica, os

programas de cálculo de engenharia sempre são comparados com casos simples para que

possam ser validados. Se o cálculo numérico submetido às mesmas condições de contorno

que um cálculo analítico apresentar resultados totalmente diversos a este, então, nem ao

menos se pode confiar no método numérico. Portanto, comparar resultados numéricos com

os analíticos é um passo importante na validação de qualquer método numérico.Tendo isso

em vista, esse capítulo se dedica a apresentar os principais elementos de sistemas de

aterramento elétrico que possuem soluções analíticas, contudo, sem demonstrações formais.

Em geral, os sistemas de aterramentos são compostos por fios ou cabos condutores

dispostos horizontalmente ou verticalmente no solo. O comprimento, o diâmetro e a

disposição destes cabos no solo são parâmetros determinantes na dispersão das correntes e,

conseqüentemente, na resistência do sistema a terra e na distribuição de potencial ao longo

do solo. Em aterramentos de maiores dimensões, vários cabos são ligados entre si de forma

a construir uma malha de aterramento, uma geometria de cálculo significativamente mais

complexa. Os casos mais simples a serem abordados são:

Haste simples disposta verticalmente

Condutor enterrado horizontalmente

Dois condutores em ângulo reto dispostos horizontalmente

Dois condutores em cruz dispostos horizontalmente no solo

As soluções analíticas são deduzidas através da aplicação das leis de Maxwell e toma-se

vantagem da simetria dos casos analisados. Assume-se que as geometrias são compostas de

condutores cilíndricos, cujos comprimentos são significativamente maiores que seus

diâmetros. Além disso, os modelos analíticos supõem que através do sistema de

aterramentos é vista uma resistividade aparente do solo, como se o solo fosse homogêneo.

Esta simplificação é razoável quando o sistema de aterramento está todo contido dentro de

uma única camada de resistividade constante. Nestes modelos, pode-se levar em conta a

13

estratificação do solo em camadas corrigindo-se o valor da resistividade aparente para uma

profundidade de interesse. Em outras palavras, a resistividade aparente é unicamente um

artifício matemático para modelar o solo em questão1. Essas assunções são válidas, pois

grande parte da energia é dissipada muito proximamente ao sistema de aterramento, onde a

resistividade pode ser considerada constante.

3.1 RESISTÊNCIA DE HASTES VERTICAIS

Uma haste delgada enterrada verticalmente no solo é o sistema de aterramento mais

simples que existe. Este tipo de sistema é largamente utilizado em residências e em postes

para aterrar o neutro do secundário dos transformadores de distribuição.

Ilustração 3-1 – distribuição das equipotenciais em uma haste simples

1 Existem diversas formas para se encontrar a resistividade aparente do solo. Cita-se, por exemplo, a fórmula de Hummel (Ref. [1]) para um solo estratificado em duas camadas.

14

A corrente injetada no topo da haste desce, se espalha pela superfície do condutor e se

propaga através do solo, formando superfícies equipotenciais elípticas que ao se

distanciarem da haste tornam-se cada vez mais esféricas. Nota-se que o campo elétrico

decresce muito rapidamente ao se distanciar da haste, isso reitera a afirmação de que grande

parte da energia dissipada se dá muito próximo da haste. As referências 1 e 2 descrevem a

equação para a resistência de uma única haste como:

−

= 14

2 r

Ln

LR a

h lπρ

Equação 3-1

Onde:

aρ é a resistividade aparente do solo [Ω.m]

L é o comprimento da haste [m]

r é o raio da haste [m]

A partir da Equação 3-1, nota-se que conforme se aumenta o comprimento da haste,

diminui-se sua resistência de aterramento. Porém, chega-se a um ponto que aumentar o

comprimento da haste já não é mais vantajoso, pois mesmo aumentando-se bastante o

comprimento da haste, a redução da resistência será muito pequena. Portanto há um

tamanho de haste economicamente ótimo para a maioria dos casos. Em geral, as hastes

comerciais têm entre 2,4 e 3,0 m de comprimento e de ½ a ¾ de polegada de diâmetro.

O caso da haste delgada também possui solução analítica para a distribuição de

potencial no solo, conforme mostra a equação 3-2.

( )d

LdLn

L

IrV a ++=

22

2l

πρ

Equação 3-2

Onde:

I é a corrente de defeito a fluir pela haste [A]

aρ é a resistividade aparente do solo [Ω.m]

L é o comprimento da haste [m]

d é a distância da haste [m]

15

3.2 RESISTÊNCIA DE CONDUTORES DISPOSTOS HORIZONTALMENTE

Várias hastes verticais podem ser interligadas de forma a reduzir a resistência do

sistema de aterramento. Em geral, essa ligação entre hastes é feita através de condutores

dispostos horizontalmente no solo a uma dada profundidade. Estes condutores de fato

contribuem ainda mais para a redução da resistência do sistema.

A interação mútua entre os condutores horizontais e os verticais é extremamente

complexa e não há solução analítica que contemple esse caso até o presente momento.

Porém, ainda assim, é interessante estudar o caso de condutores horizontais enterrados para

os fins de validação dos programas e para alguns casos práticos como a resistência de cabos

contra-peso de torres de transmissão.

Ilustração 3-2 – condutor disposto horizontalmente

16

Em [1], apresenta-se a seguinte equação para um condutor singelo enterrado no solo:

−+−

= Ll

L

p

rp

Ln

LR a 22

2

2

πρ

Equação 3-3

Onde:

aρ é a resistividade aparente do solo [Ω.m]

L é o comprimento do condutor [m]

r é o raio do condutor [m]

p é a profundidade [m]

Há também outras formas de se dispor condutores no solo. Por exemplo, dois

condutores formando um ângulo reto, ou então uma cruz. Como nas figuras abaixo. A

referência [1] propõe solução analítica para estes casos (Equação 3-4 e Equação 3-5).

Ilustração 3-3 condutores dispostos em L e em cruz

17

Equação para dois condutores enterrados horizontalmente em ângulo reto:

++−

= Ll

L

p

rp

Ln

LR a 4292,02373,0

4

2

πρ

Equação 3-4

Equação para dois condutores enterrados horizontalmente em cruz:

+−+

= Ll

L

p

rp

Ln

LR a 284,4526,1

4

2

πρ

Equação 3-5

A referência [1] ainda mostra que, dado um comprimento de cabo para ser enterrado

horizontalmente, a configuração de menor resistência é a do condutor singelo, seguido pelo

ângulo reto e depois pela cruz. Isso se explica em razão da interferência e da redução do

comprimento efetivo das pontas. De acordo com o exemplo da referência, dispõe-se de 60

m de um condutor com 6 mm de diâmetro, a ser enterrado a 60 cm de profundidade em um

solo com resistividade aparente de 1000 Ω.m. A Tabela 3-1 mostra os resultados.

Tabela 3-1 – comparação entre geometrias Configuração Resistência [Ω.m]

Um único condutor 35,0

Dois condutores em ângulo reto 64,8

Dois condutores em cruz 73,2

Estes quatro sistemas de aterramento que possuem solução analítica serão utilizados

como base para a validação do programa de elementos finitos no final deste trabalho.

18

3.3 RESISTÊNCIA DE UMA MALHA DE ATERRAMENTO SEGUNDO A

IEEE STD 80

Em sistemas de potência, exige-se uma baixíssima resistência de terra, tanto por

critérios técnicos, relacionados à proteção do sistema, como também por critérios de

segurança a vida. Os sistemas de aterramento devem ser dimensionados de tal maneira que

a distribuição de potencial elétrico causado pela falta não torne um operador próximo

sujeito a tensões elevadas. Em outras palavras, a diferença de potencial entre suas

extremidades como mãos e pés, etc, não pode ultrapassar um valor crítico, acima do qual

causaria fibrilação do coração.

Para satisfazer as duas condições, a baixa resistência a terra e baixos potenciais de toque

e de passo, as subestações elétricas costumam utilizar uma malha de terra composta de

cabos horizontais em quadriculados e hastes verticais conectadas a alguns (ou todos) os nós

da malha. A Ilustração 3-4 ilustra uma típica malha de terra e seus parâmetros geométricos.

As malhas de terra são muito eficientes e econômicas, e são utilizadas em quase todas as

subestações.

Ilustração 3-4 – parâmetros geométricos de uma malha retangular.

19

Para a malha de terra não há solução analítica. Mesmo que se compusesse a malha

como a resistência paralela de vários condutores enterrados horizontalmente e hastes

verticais, o resultado seria significativamente menor do que o real em razão do desprezo

das intereferências. Logo, não há modelo analítico que possa representá-la. O que se faz

desde há muito tempo é utilizar uma fórmula semi-empírica, descrita na norma IEEE Std 80

[2]. De fato, as malhas das subestações utilizadas como exemplo neste trabalho foram

calculadas de acordo com esta norma. A norma cita a fórmula de Sverak para uma primeira

estimativa da resistência da malha.

Fórmula de Sverak:

+++=

mmmtam

AhALR

/201

11

20

11ρ Equação 3-6

Onde:

Rm é a resistência da malha

ρa é a resistência aparente do solo

Lt é o comprimento total de cabos e hastes

Am é a área ocupada pela malha

hm é a profundidade da malha

Neste capítulo, demonstraram-se as soluções analíticas para alguns sistemas simples de

aterramento feitos por condutores e hastes e o método proposto pela norma IEEE Std 80

para malhas de subestações. Várias páginas deste trabalho foram dedicadas a estas soluções

analíticas e semi-empíricas para que depois os resultados numéricos do programa possam

ser validados.

Aqui se encerram as discussões sobre métodos analíticos, semi-empíricos,

estratificação, medida e modelagem do solo, para dar lugar ao cerne deste trabalho: a

utilização dos métodos numéricos para o cálculo de sistemas de aterramento elétricos.

20

4. MÉTODO DOS ELEMENTOS FINITOS

4.1 MÉTODOS NUMÉRICOS PARA ELETROMAGNETISMO

Os problemas de aterramento são típicos problemas de eletromagnetismo, que é uma

ciência regida por equações diferenciais vetoriais, e eles têm sido resolvidos

numericamente por três tipos de métodos [3]: os baseados em teoria de circuitos e linhas de

transmissão, aqueles baseados no cálculo de campos eletromagnéticos e os híbridos.

O método que se baseia na teoria de circuitos é o Método das Linhas de Transmissão ou

TLM (Transmission Line Method), que se utiliza da analogia existente entre as equações

que regem os campos eletromagnéticos e as equações de propagação de ondas em uma

linha de transmissão. Este método cria correspondência física ao substituir o domínio de

estudo por um reticulado de linhas de transmissão onde as correntes e tensões nos seus nós

representam os campos magnético e elétrico, respectivamente. As linhas de transmissão

podem ser a parâmetros concentrados ou distribuídos e o sistema é resolvido da mesma

forma que os sistemas elétricos de potência, através de matrizes de rede.

O TLM se baseia, também, no modelo de propagação de onda de Huygens. Dada uma

perturbação em um ponto qualquer do domínio de estudo, leva um tempo até que os outros

pontos tomem ciência. Isto aproxima ainda mais este modelo da realidade física. O TLM é

vantajoso ao lidar com situações que envolvam fenômenos de propagação (alta freqüência),

como uma descarga atmosférica, por exemplo. Porém, em baixas freqüências, o método

tende a ficar impreciso. O TLM tem algumas desvantagens, por exemplo, sua implantação

não é trivial. Além disso, é muito difícil de modelar meios anisotrópicos, o que torna muito

trabalhosa a implantação de PMLs. Este método foi utilizado raramente para simulação de

sistemas de aterramento em regime estacionário.

21

Os métodos com base na teoria eletromagnética são mais robustos que aqueles baseados

na teoria de circuitos, pois admitem menos hipóteses aproximadoras [3] e também há

também o fato de que são mais adequados a problemas em baixa frequência. Ainda assim,

há muitos métodos baseados na teoria eletromagnética e este trabalho discorrerá muito

brevemente sobre eles e suas aplicações.

Podem-se dividir os métodos numéricos para resolução de problemas de equações

diferenciais parciais, ou problemas de contorno, em duas grandes classes: métodos

integrais, ou de fronteira, e os métodos diferenciais, ou de domínio.

Os métodos de fronteira são voltados para a solução de equações integrais e tratam de

encontrar condições de contorno para estas equações que se adequem às condições de

contorno na fronteira. Numa etapa posterior, de pós-processamento, a distribuição interna

dos campos é determinada analiticamente através das equações integrais. Essa característica

peculiar dos métodos de fronteira faz com que eles sejam capazes de lidar com geometrias

infinitas ou semi-infinitas e, por essa razão, têm sido largamente utilizados em problemas

de antenas. O Método dos Momentos e o Método dos Elementos de Fronteira (ou

Contorno) são exemplos clássicos de métodos integrais.

Há dois grandes inconvenientes nestes métodos de fronteira que limitam a sua

aplicação. Um deles é a limitação a problemas cujos meios são homogêneos, isotrópicos e

lineares e o outro é que esses métodos levam a matrizes cheias. Essas matrizes crescem

com o quadrado do tamanho do problema e passam a exigir grandes quantidades de

memória para problemas relativamente pequenos. Já os métodos de domínio levam a

matrizes esparsas, normalmente de bandas, que crescem quase que linearmente com o

tamanho do problema. Por essa razão os métodos de domínio são geralmente muito mais

eficientes que os de fronteira para problemas relativamente grandes e, portanto, têm mais

aplicações a estudos eletromagnéticos, em especial ao estudo do aterramento.

Os métodos diferenciais partem da discretização do domínio de estudo em uma malha

para então discretizar as equações diferenciais que regem o problema, de forma que elas

possam ser escritas como equações algébricas lineares, passíveis de serem resolvidas por

um computador. Essas equações lineares podem ser solucionadas iterativamente, como

22

pelos métodos de Gauss-Seidel, ou então deterministicamente como pelo método da

Eliminação de Gauss com pivotamento, etc.

Dos diversos métodos de domínio aplicáveis a problemas eletromagnéticos, dois são de

grande interesse para o estudo de aterramentos: o Método das Diferenças Finitas (MDF) e o

Método dos Elementos Finitos (MEF).

O Método das Diferenças Finitas divide o domínio a ser estudado em uma malha

estruturada de nós, aproxima a equação diferencial que rege o fenômeno e as condições de

contorno em equações lineares algébricas mais simples (chamadas equações de diferenças

finitas) e então resolve esse sistema algébrico de equações aproximadas. Se o número de

nós da malha for suficientemente grande, isto é, se a distância entre os nós for muito

pequena, no limite a solução das equações algébricas tenderá a ser a mesma das equações

diferenciais.

Quando se trata de métodos numéricos para eletromagnetismo, o MDF é geralmente

utilizado para se resolver problemas no domínio do tempo, o tão chamado FDTD (Finite

Difference Time Domain). É um método direto, de formulação muito simples e fácil de ser

implantado. As equações de Maxwell são discretizadas a partir de aproximações de

diferenças centrais em relação às derivadas no tempo e no espaço. Logo, essas equações

algébricas de diferenças finitas (não mais equações parciais diferenciais), são resolvidas

progressivamente para cada instante de tempo, partindo-se de um instante zero e condições

iniciais conhecidas. Por exemplo, para um dado instante de tempo, num primeiro momento

resolve-se o campo elétrico ao longo do domínio de estudo e num segundo momento

resolve-se então o campo magnético. Então, parte-se para o próximo instante de tempo.

Repetem-se os cálculos até que o período a ser analisado se complete ou então até que o

sistema atinja o regime permanente.

Um dos pontos mais fortes do FDTD é que ele pode abranger uma gama de freqüências

muito grande, desde freqüências muito baixas da ordem de alguns Hertz a freqüências

elevadíssimas da ordem dos GHz. Esse método tem sido utilizado com sucesso na

resolução de problemas de aterramento onde se está interessado no transitório, como a

condução de uma descarga atmosférica para a terra. Porém, quando freqüências muito altas

23

estão presentes, o espaçamento entre os nós, tanto no espaço quanto no tempo, deve ser

consideravelmente pequeno, para que pequenos comprimentos de onda possam ser

devidamente representados. Isso cria malhas de tamanhos exorbitantes, muito pesadas

computacionalmente e que pode levar a tempos proibitivos de simulação.

Outra grande desvantagem do MDF é que para facilitar a aproximação por diferenças

finitas, normalmente adota-se uma malha espacial reticulada (estruturada), feita através

cubos. Este tipo de malha não é conveniente para algumas geometrias curvas, que são

muito comuns, como cilindros (modelo de um condutor, por exemplo). Esse é o fato

gerador de muitos erros do método, exigindo uma discretização consideravelmente mais

fina para tentar acompanhar os contornos mais curvos. Apesar disso, o MDF em especial o

FDTD, tem sido utilizado para simulação de transitórios eletromagnéticos em sistemas de

aterramento.

O método dos elementos finitos (MEF) baseia-se na discretização do domínio de estudo

em diversos subdomínios menores (ou elementos), onde os parâmetros a serem analisados

em seu interior, como tensões, temperaturas, potenciais, etc., são aproximados por alguma

função interpoladora a partir dos valores em seus nós. Há formulações mais avançadas do

MEF que utiliza elementos de arestas, ao invés de elementos de nós, que têm sido

empregados com muito sucesso por [3] em problemas de altas frequencias que envolvem

propagação. Neste trabalho, limitar-se-á apenas ao MEF nodal.

Ainda que o MEF seja mais complexo, menos intuitivo, conceitualmente mais

elaborado e significativamente mais difícil de ser implantado que o MDF, ele é

consideravelmente mais versátil que este, pois os elementos finitos não precisam ser

necessariamente iguais. Isto faz com que domínios de estudo com geometrias complexas

possam ser analisados com maior precisão. Além disso, podem-se estabelecer mais

elementos em regiões de interesse, onde os campos variam rapidamente ou o raio de

curvatura é pequeno, e diminuir o número de elementos em regiões distantes.

24

4.2 APLICAÇÃO DO MEF EM SISTEMAS DE ATERRAMENTO

4.2.1 Formulação Matemática

Nos jargões do eletromagnetismo, o estudo do aterramento elétrico em regime

permanente se caracteriza como um problema de eletrocinética. De outra maneira, trata-se

de um problema em que cargas elétricas em movimento constante geram campos elétricos e

magnéticos invariantes no tempo.

As equações que regem a eletrocinética podem ser facilmente derivadas a partir das

equações de Maxwell ao se levar suas derivadas parciais com relação ao tempo a zero e

assumir que a carga líquida em um ponto é nula. A tabela X abaixo resume as equações de

Maxwell em sua forma geral e para o caso particular da eletrocinética.

Tabela 4-1– Equações para o eletromagnetismo geral e eletrocinética nas suas formas diferenciais

Forma Geral Eletrocinética

Lei de Faraday t

BE

∂∂−=×∇r

r 0=×∇ E

r

Lei de Ampère t

DJH

∂∂+=×∇r

rr JH

rr=×∇

Lei de Gauss ρ=⋅∇ Dr

0=⋅∇ Dr

Lei de Gauss para o Magnetismo 0=⋅∇ Br

0=⋅∇ Br

Lei da Continuidade t

J∂∂−=⋅∇ ρr

0=⋅∇ Jr

Relação Constitutiva EJrr

σ= EJrr

σ=

25

A lei da continuidade assume a indestrutibilidade da carga, pois se o fluxo líquido de

cargas em um dado ponto do espaço não for nulo, então este ponto estará acumulando

cargas elétricas. Porém, no regime eletrocinético, a quantidade de carga em um ponto não

varia com o tempo, em outras palavras, não há acumulo de cargas. Por sua vez, a equação

constitutiva relaciona uma propriedade física do material a variáveis elétricas.

A partir da lei da continuidade e da relação constitutiva deduz-se a equação do

problema deste estudo. Para tanto, relaciona-se o campo elétrico com o potencial elétrico V:

( )VEJ ∇−== σσrr

( ) 02 =∇−=∇−⋅∇=⋅∇ VVJ σσr

Destas relações, fica claro que a equação que rege o aterramento elétrico em regime

permanente é a equação conservativa de Laplace. Há também outras condições de fronteira

que devem ser obedecidas para garantir a unicidade à solução. A superfície fechada

ND Γ+Γ=Γ que delimita o domínio de estudo Ω é dividida em duas fronteiras, onde se

aplicam diferentes condições de contorno, como mostrado na Ilustração 4-1.

ΓN

ΓD

ΓD

Ilustração 4-1 domínio de estudo truncado

26

A fronteira DΓ , que foi truncada, deve representar o ponto de terra remoto (infinito) e,

portanto, nela se impõe potencial elétrico nulo. Em tese, o potencial elétrico a uma

distância finita da fonte de campo tem um valor finito não nulo, porém é razoável assumir

que o potencial 0V neste ponto de fronteira seja praticamente zero desde que ele esteja a

uma distância razoável da fonte de campo. A essa imposição denomina-se condição de

Dirichlet.

Quanto à fronteira NΓ , que representa a superfície da subestação, deve-se impor a

condição de von Neumann, pois se assume que não há corrente elétrica atravessando-a. É

uma aproximação muito conveniente e adequada, pois se considera que o ar é um isolante

perfeito. Essa hipótese também é muito razoável, já que a condutividade do ar é

significativamente menor do que a do solo comum.

Em suma, as equações X, Y e Z resumem a formulação matemática que governa o

aterramento elétrico em regime estacionário.

02 =∇ Vσ em Ω Equação 4-1

0VV = em DΓ Equação 4-2

0ˆ =⋅ nJr

, ou 0ˆ =⋅∇ nV , ou ainda 0ˆ

=∂∂

n

V em NΓ Equação 4-3

Onde n é um vetor normal à fronteira de superfície.

Então o método dos elementos finitos é usado para se resolver estas equações

simultaneamente e este processo se inicia pela divisão do domínio de estudo.

27

4.2.2 Divisão do Domínio de Estudo

O primeiro passo na implantação do MEF para se resolver um problema consiste na

divisão do domínio de estudo em regiões menores, em elementos, de tal forma que eles

sejam conectados em seus vértices (nós). O MEF calculará o potencial elétrico nos nós da

malha criada, e o aproximará dentro de um elemento através de uma função que será

discutida mais adiante. O importante neste primeiro momento é que a malha seja adequada

ao tipo de problema a ser resolvido. Isso significa que deve haver mais elementos em

regiões com maior variação (gradiente) de potencial e elementos menores nas regiões de

pequeno raio de curvatura para se conformarem melhor ao perímetro, de forma melhorar a

precisão numérica. No caso do estudo de aterramento, isso significa que deve haver muitos

elementos no entorno dos condutores, tanto pelo fato geométrico quanto pela magnitude do

campo elétrico.

O processo de divisão do domínio de estudo pode ser feito manualmente, o que é muito

trabalhoso para geometrias muito grandes e complexas, ou, como é mais comum, através de

um gerador de malha que divide automaticamente o domínio a ser estudado. Neste trabalho,

os autores utilizaram o GMSH [7], um software livre, para criar as malhas, que será citado

mais adiante.

Há diversos tipos de malhas, como as malhas estruturadas, que possuem regularidade na

formação de seus elementos e as malhas não estruturadas, muito irregulares, mas que se

conformam melhor em geometrias mais complexas. Podem-se utilizar os mais diversos

elementos para se discretizar o domínio, cada qual com suas características. O tipo de

elemento influi no tempo de processamento do problema e na precisão dos resultados. A

Ilustração 4-2 mostra alguns elementos.

28

Elementos bidimensionais

Elementos tridimensionais

Elemento unidimensional

Ilustração 4-2 elementos típicos para o método dos elementos finitos

O GROUND-3D utiliza elementos de uma, duas e três dimensões. Os elementos

unidimensionais são usados para representar os condutores da malha, já que seus

comprimentos são muito mais significativos que seus diâmetros. Os elementos

bidimensionais são triângulos usados para representar camadas cuja espessura são

insignificantes frente ao tamanho do domínio de estudo, como a camada superficial de brita

da subestação. Por fim, os elementos tridimensionais são prismas de base triangular usados

para representar o solo em si.

29

4.2.3 Funções Interpoladoras e Funções de Forma

Depois de se ter definido a malha, define-se uma equação que aproxime o potencial

elétrico ao longo de um elemento. Essas equações são tipicamente funções polinomiais, que

são mais simples de derivar. Essas funções interpolam o valor do potencial elétrico dentro

do elemento a partir dos potenciais em seus nós. Dessa forma, ao se somar as funções

interpoladoras de todos os elementos do domínio têm-se uma aproximação do parâmetro de

estudo (potencial elétrico) ao longo de todo o domínio (solo).

Então, se N for o número total de elementos no domínio, ( )zyxVe ,, for uma função

interpoladora do potencial elétrico para o elemento e, tendo um valor qualquer dentro desse

elemento e valendo zero fora dele, então se pode escrever que o potencial elétrico em

qualquer ponto do domínio ( )zyxV ,, como:

( ) ( )∑=

≅N

ee zyxVzyxV

1

,,,, Equação 4-4

A função interpoladora ( )zyxVe ,, pode assumir diversas formas. Ela define em boa

parte a precisão da simulação e para cada tipo de elemento há uma função interpoladora

adequada. Por exemplo, para um elemento tetraédrico que possui apenas 4 nós (em seus

vértices), uma função interpoladora típica é:

( ) dzcybxayxVe +++=, Equação 4-5

A Equação 4-5 é chamada de função interpoladora de primeira ordem e é utilizada em

elementos de primeira ordem, porque os parâmetros x, y e z são elevados à unidade. O

triângulo de três nós e o quadrilátero de quatro nós da Ilustração 4-2 são exemplos de

elementos de primeira ordem.

Em algumas aplicações, onde se exige maior precisão dos cálculos, elementos de

segunda ordem são utilizados, como o triângulo de seis nós e o quadrilátero de oito nós da

Ilustração 4-2. Nestes elementos são utilizadas funções interpoladoras que possuem termos

30

quadráticos de forma a possibilitar uma melhor interpolação. Porém, em razão do maior

número de nós na malha, o tempo de simulação aumenta consideravelmente, a despeito da

maior precisão obtida.

Para os elementos lineares, que possuem uma função interpoladora como a da Equação

4-5, o gradiente de potencial ao longo do elemento é constante, isto é, o campo elétrico ao

longo do elemento é constante:

zdycxbV ˆˆˆ ++=∇ Equação 4-6

zEyExEE zyx ˆˆˆ ++=r

Equação 4-7

Os coeficientes de interpolação para um elemento qualquer são determinados em função

dos potenciais elétricos nos nós deste mesmo elemento. Por um artifício matemático, pode-

se escrever a interpolação no interior do elemento como um somatório dos potenciais

elétricos nodais ponderados por uma função ( )zyxi ,,α . Mais adiante, essas funções serão

úteis para simplificar o problema de se unir os elementos em uma única equação global.

( ) ( )∑=

=n

ieiie VzyxyxV

1

,,, α Equação 4-8

A esses fatores de ponderação denominam-se funções de forma do elemento (element

shape function) e eles nada mais são que funções de interpolação de ordem igual à do

elemento. Essas funções de forma têm algumas características especiais, que reforçam a sua

natureza ponderadora: elas assumem a unidade quando o ponto de interesse (x,y,z) coincide

com o nó a que ela se refere e assume zero nos outros nós.

31

4.2.4 Formulação Matricial do Problema

O método de Rayleigh-Ritz, que se baseia no cálculo variacional, por esta razão também

é conhecido como método variacional, resolve um problema de contorno ao aproximar a

solução com uma combinação linear finita de funções mais simples. O cálculo variacional

tem por objetivo otimizar (normalmente minimizar) uma classe especial de funções

denomidas funcionais.

A equação de Laplace, que governa o problema do aterramento elétrico no regime

estacionário, é conservativa e sua solução representa um estado de mínima energia para o

sistema. O problema então pode ser resolvido se houver uma forma algébrica de se

representar a energia do sistema em função dos potenciais nodais e, então, determinar a

combinação de potenciais que minimiza essa energia, mas que também satisfaça as

condições de contorno.

Um funcional para a equação de Laplace que rege o problema é aquele que descreve a

energia como aquela dissipada por efeito Joule no solo devido à circulação de corrente.

Partindo-se da interpretação física desse funcional, escreve-se a densidade volumétrica de

potência dissipada como função da condutividade do solo e do quadrado da magnitude do

campo elétrico local.

2

eeeeee EEEEJprrrrr

σσ =⋅=⋅= Equação 4-9

Integrando-se a densidade volumétrica de potência dissipada por efeito Joule ao longo

de todo o elemento, encontra-se a potência dissipada pelo próprio elemento:

∫∫∫∫ ∇== dSVdSEP eee

22σσ

r Equação 4-10

Calculando-se o gradiente do potencial V, obtém-se o campo elétrico para qualquer

ponto de interesse e, dessa forma, a potência dissipada por efeito Joule. Partindo-se das

aproximações para o potencial elétrico no domínio de estudo, vem:

32

( ) ( )∑∑==

∇=

∇=∇=3

1

3

1

,,,,i

ieii

ieiee zyxVzyxVVE ααr

Equação 4-11

Substituindo-se a Equação 4-11 na Equação 4-10, encontra-se uma forma de se

expressar a potência dissipada no interior do elemento como uma função dos potenciais

nodais.

[ ] eji j

jieie VdSVP ∑∑ ∫= =

∇⋅∇=3

1

3

1

αασ Equação 4-12

Definindo-se o termo entre colchetes como:

( )∫ ∇⋅∇= dSG ji

eij αασ Equação 4-13

Pode-se representar a potência dissipada por efeito Joule no interior de um elemento na

forma matricial, como função das tensões nodais:

[ ] ( )[ ][ ]eet

ee VGVP = Equação 4-14

Onde t denota a transposição do vetor de potenciais elétricos e ( )[ ]eG é a matriz dos

coeficientes do elemento, também conhecida como matriz de rigidez, devido à aplicação

pioneira do MEF em análises estruturais. A ordem (tamanho) dessa matriz depende do tipo

de elemento usado e sua ordem é igual ao número de nós do elemento.

Os gradientes das funções de forma, por serem polinômios, são fáceis de calcular e

implantar computacionalmente. Logo, tem-se a potência dissipada em um elemento em

função dos potenciais nodais Ve. O próximo passo é montar os elementos finitos num

processo conhecido como assemblagem (do inglês assembly), relacionando as equações de

todos os elementos em um único sistema global. Por conta disso, a matriz desse sistema é

denominada matriz global, e é dela que se trata a próxima seção.

33

4.2.5 A Matriz Global e a Potência Dissipada

Assim que os coeficientes das matrizes de todos os elementos forem calculados,

prossegue-se para a montagem da matriz global. A partir deste ponto, todas as numerações

dos nós serão globais e não mais locais, como havia sido feito até então. A matriz global

[ ]G é obtida a partir do fato que a distribuição de potencial elétrico deve ser contínua entre

as fronteiras dos elementos. Demonstra-se que um coeficiente qualquer ijG é soma das

contribuições de ( )epqG de cada elemento onde os nós locais p e q são os mesmos que os nós

globais i e j.

A matriz global também é conhecida como matriz de condutâncias, pois a unidade de

seus coeficientes é de condutância (Siemens). Existe uma interpretação física da matriz

global, ela pode ser enxergada como a matriz de rede de um sistema puramente resistivo,

onde os nós dos elementos são os nós dos circuitos elétricos e as arestas dos elementos são

condutâncias que ligam os nós.

Com a matriz global montada, calcula-se a potência dissipada pelo sistema:

[ ] [ ][ ]VGVPP tN

ee ==∑

=1

Equação 4-15

Onde [V] é o vetor de potenciais nodais com referência global:

[ ]

=

nV

V

V

V

V

M

3

2

1

Equação 4-16

Há algumas características muito interessantes que a matriz global possui que acabam

tornando o MEF apelativo. Primeiramente, a matriz global é simétrica, isto é, jiij GG = . A

simetria traz algumas vantagens em operações matriciais, como o cálculo de determinantes.

34

Em segundo lugar, a matriz global é largamente esparsa para sistemas grandes e pode se

tornar uma matriz de bandas se os nós forem cuidadosamente numerados. Essas

características de esparsidade e banda favorecem enormemente os cálculos e reduz

significativamente o tempo de computação. Por fim, a matriz é singular, de modo que é

impossível invertê-la numericamente, tornando o problema de natureza iterativa.

4.2.6 A Solução do Sistema

O sistema é resolvido utilizando-se o método dos gradientes conjugados, um método

iterativo que resolve o sistema através da minimização de um funcional quadrático. Para

tratar do mau condicionamento do sistema utiliza-se uma técnica de pré-condicionamento

denominada Decomposição Incompleta por Cholesky. A associação desses dois métodos

forma o ICCG (Incomplete Cholesky Conjugate Gradientes) ou Método dos Gradientes

Conjugados pré-condicionado por Decomposição Incompleta da Cholesky.

O método dos gradientes conjugados é muito eficiente para sistemas com matrizes

simétricas com termos positivos. Basicamente, o método opera gerando seqüências de

vetores de soluções aproximadas, vetores de resíduos e vetores de direções de procura, que

são usados para atualizar os próximos vetores solução aproximada e resíduos. Esse método

é muito interessante do ponto de vista dos requisitos computacionais, pois exige o

armazenamento na memória de apenas três vetores, além dos termos não nulos da matriz

global.

O pré-condicionamento de Cholesky trata-se de uma transformação matricial que

aproxima a matriz do sistema à matriz identidade, que continua altamente esparsa,

acelerando a convergência do método dos gradientes conjugados.

O ICCG é um método muito eficiente e rápido, que se mostrou muito adequado para a

solução de problemas de aterramento.

35

4.2.7 Resistência do Sistema de Aterramento

Após a solução do problema pelo ICCG, conhece-se o potencial elétrico em cada nó da

malha de elementos finitos. Substituindo-se o vetor potencial elétrico V na equação 4-15

obtém-se a potência dissipada no solo por efeito Joule. A partir da elevação de potencial no

ponto de defeito Vd e da potência dissipada no solo é possível avaliar a resistência aparente

do aterramento.

P

VR d

2

= Equação 4-17

Na realidade, as versões mais recentes do GROUND-3D, como a utilizada pelos autores

para a realização deste trabalho, não se utilizam mais da alimentação do sistema em tensão

e sim em corrente, logo:

2I

PR = Equação 4-18

Com a solução do problema (distribuição de potencial elétrico) também é possível

prever o potencial de passo máximo a que estaria sujeito um operador dentro da subestação

quando da ocorrência de uma falta. Este trabalho não se preocupou em analisar o potencial

de passo, pois os autores julgaram que a resistência a terra do sistema já é um parâmetro

representativo da precisão da simulação.

36

4.3 PML (PERFECTLY MATCHED LAYERS)

No início deste capítulo foi citado que há um erro inerente aos métodos numéricos

devido à discretização do domínio. Porém, este não é o único erro inerente ao método

numérico, pelo menos quando o problema a ser resolvido trata de eletromagnetismo.

Uma carga elétrica produz um campo elétrico e, se em movimento, um campo

magnético que se estendem ad infinitum. Em tese, o domínio de estudo deveria ser infinito

para contemplar esse fato. Ainda que seja impossível representar um domínio de estudo

infinito computacionalmente pelo método dos elementos finitos, pode-se truncar o domínio

a uma distância razoável da região de interesse e assumir que o erro devido ao truncamento

é muito pequeno e aceitável.

Os campos elétricos e magnéticos a uma grande distância de suas fontes decaem com o

inverso da distância porque podem as fontes podem ser vistas como uma única fonte

pontual. Logo, é razoável truncar o domínio de estudo se ele for dezenas a centenas de

vezes maiores que a região de interesse, pois os campos elétricos e magnéticos nas suas

extremidades serão muito próximos de zero.

A princípio, pode parecer que mais exata será a solução numérica ao se aumentar o

domínio de estudo. Porém, a memória computacional é limitada e o aumento do domínio

exigirá que a discretização seja menor, isto é, que os elementos finitos sejam maiores.

O truncamento traz consigo outro problema ainda mais grave, denominado reflexão

numérica, que surge devido à descontinuidade na interface, isto é, na fronteira artificial.

Uma maneira de se contornar isso é adequar as condições de contorno na fronteira ou então

utilizar o método dos elementos finitos com formulação significativamente diferente e

específica. Por exemplo, há uma formulação híbrida conhecida como MEF-MEC (Método

dos Elementos de Contorno), porém esta reduz substancialmente a esparsidade da matriz

global, que é uma das grandes vantagens do MEF original [3].

37

Os problemas oriundos do truncamento, tanto o elevado tamanho do domínio do estudo

como a reflexão numérica, podem ser significativamente atenuados com uma solução

proposta em 1994 por Bérenger [9], denominada PML – Perfectly Matched Layers, ou

camadas perfeitamente casadas.

O PML foi derivado a partir das equações de Maxwell modificadas com a aplicação de

uma transformação de coordenadas que “esticam” os eixos cartesianos. Essa abordagem,

demonstrada em [8], é muito pouco intuitiva. Há outra abordagem, mais simples e com uma

correspondência física.

O PML pode ser considerado como uma camada ou mais de material anisotrópico de

altas perdas que absorve completamente qualquer campo incidente na interface entre o

PML e o domínio de estudo, independentemente de sua freqüência, magnitude, polarização

ou ângulo de incidência.

Ilustração 4-3 truncamento do domínio por material absorvedor anisotrópico PML

Sendo o PML anisotrópico, as propriedades do material, como a condutividade elétrica,

σ , dependem do sentido do vetor densidade de corrente, mais especificamente, da posição

da camada PML em relação à região de interesse. Ao se escrever as equações de Maxwell

em sua forma diferencial se inclui a anisotropicidade do material através de um tensor para

a condutividade elétrica, σ , na equação da continuidade.

38

EJ ⋅= σ Equação 4-19

Onde o tensor de condutividade é escrito como:

=

zz

yy

xx

σσ

σσ

00

00

00

Equação 4-20

Como já mencionado, o PML foi deduzido através de uma transformação de

coordenadas que esticavam ou contraiam os eixos cartesianos dentro da camada

absorvedora e que encontrou uma correspondência física em um material anisotrópico.

Então, o tensor da condutividade elétrica pode ser escrito em função de um único tensor

Λ que contempla a transformação de coordenadas.

Λ= σσ Equação 4-21

Onde σ representa a condutividade elétrica do material da fronteira que está em contato

direto com o PML. Por fim, o tensor Λ é escrito em função de “esticadores de

coordenadas”, xxs , yys e zzs , para o caso cartesiano.

=Λ

zz

yy

xx

s

s

s

00

00

00

Equação 4-22

Esses parâmetros “esticadores de coordenadas”, por sua vez, dependem do plano em

que o PML está localizado, em outras palavras, dependem da direção normal deste plano

que atenuará com máxima intensidade as correntes incidentes. No caso das simulações

realizadas neste trabalho, as camadas de PML são sempre ortogonais às direções

cartesianas x , y e z . Por exemplo, se a camada de PML for ortogonal ao eixo z, a

transformação de coordenadas poderia ser feita através das relações [3]:

ass yyxx == Equação 4-23

aszz1=

Equação 4-24

39

No limite, se as dimensões dos elementos tendessem a zero, o PML funcionaria

perfeitamente. Intuitivamente, quão menores forem as dimensões dos elementos próximos à

interface do PML, melhor serão as propriedades de absorção. Isto traz alguns reveses, como

o aumento do número total de elementos. Além disso, o coeficiente de reflexão na interface

é inversamente proporcional à espessura do PML. Uma espessura muito pequena pode não

atenuar os campos suficientemente, mas uma espessura muito grande pode comprometer o

tamanho da malha. Não fez parte do escopo deste trabalho determinar as melhores

dimensões e parâmetros para o PML e preferiu-se utilizar os valores propostos em [3].

40

4.4 GMSH

Como já discutido, a geração de malha se faz necessária em um grande número de

métodos numéricos, inclusive o dos elementos finitos. As características das malhas afetam

significativamente os resultados obtidos, logo é muito importante que elas sejam elaboradas

com cuidado, com melhor discretização aonde se esperam maiores gradientes, com os tipos

e números de elementos mais adequados para o estudo em questão. Criar boas malhas é

pré-requisito para simulações representativas, e criá-las se tornou uma especialização da

engenharia. Surgiram no mercado inúmeros aplicativos de geração de malhas, cada qual

com suas capacidades e limitações.

Há quase uma década, surgiu o GMSH, um aplicativo de distribuição livre capaz de

gerar malhas tridimensionais automaticamente ou de forma assistida, com recursos de pré e

pós-processamento. Isso significa que a definição da geometria, a geração da malha e as

apresentações dos resultados podem ser feitas através deste mesmo programa. A Ilustração

4-4 e a Ilustração 4-5 mostram a utilização do GMSH para a definição da geometria de uma

máquina síncrona de pólos salientes e a Ilustração 4-6 mostra a geração sua malha.

O GMSH é constituído por quatro módulos principais: descrição geométrica, geração da

malha de elementos finitos, resolução do método numérico e pós-processamento. Todos os

comandos envolvendo estes módulos são prescritos ou com a interface gráfica para o

usuário (GUI – Graphical User Interface) ou através de arquivos de dados que utilizem a

linguagem de programação própria do GMSH. Contudo, o LMAG possui o seu próprio

programa para resolução de problemas tridimensionais de aterramento, o GROUND 3D, e

para manter a comunicação entre este e o GMSH, optou-se pela programação em arquivos

de dados. Essa escolha mostrou-se vantajosa, pois facilita as possíveis alterações nos

códigos de definição da geometria e de geração da malha.

41

Ilustração 4-4 – geração da geometria

Ilustração 4-5 – geração da malha

Ilustração 4-6 – adensamento de volumes nos pólos e no entreferro

42

Há algum tempo atrás, quando o GROUND 3D ainda não era ligado ao GMSH, a

definição da geometria e a geração da malha se dava manualmente, o que consumia muito

tempo. A utilização do GMSH possibilitou a definição da geometria por parametrização,

através de scripts, dessa forma consegue-se a geração automática dos domínios de estudos e

suas malhas de elementos finitos.

Ilustração 4-7 – distribuição de pressão nas palhetas de um estágio de uma turbina para geração termoelétrica.

Assim que a simulação é concluída pelo GROUND 3D, a interpretação dos resultados é

feita através do módulo de pós-processamento do GMSH. Ele associa a cada nó da malha

variáveis, escalares ou um vetoriais, como o potencial e campo elétrico, por exemplo.

Dessa forma, o programa é capaz de gerar gráficos de forte apelo visual que facilitam

bastante o entendimento e a interpretação dos resultados. Por exemplo, podem-se gerar

mapas de cores, mapas de campo, curvas de valor constante, linhas equipotenciais, linhas

de corrente, etc. A Ilustração 4-7 provê um exemplo das capacidades do programa.

43

Para os casos analisados neste trabalho, primeiramente geraram-se as geometrias a

serem analisadas e dividiram-se as regiões do domínio. Num próximo passo, o GMSH fez a

discretização do domínio de estudo de acordo com parâmetros especificados pelos autores.

Todas as malhas de elementos finitos possuem apenas elementos prismáticos de base

triangular para modelar o solo e unidimensionais para modelar os condutores. Isto porque o

GROUND 3D só está habilitado a trabalhar com estes tipos de elementos, pois tipos

diferentes exigiriam funções de forma diferentes.

A Ilustração 4-8 mostra, do lado esquerdo, a geração da geometria onde se pode ver os

subdomínios e, do lado direito, o domínio já discretizado em elementos prismáticos.

Ilustração 4-8 – geometria criada no Gmsh e sua malha de elementos finitos

A utilização do GMSH viabilizou este trabalho, pois dispensou os autores de se

preocuparem com a programação de um gerador de malhas e de um pós-processador para

se preocuparem tão somente com a aplicação do método numérico aos sistemas de

aterramento. Cabe notar que o próprio GMSH levou muitos anos até amadurecer. Maiores

informações sobre este programa podem ser encontradas na referência [7].

44

5. APLICAÇÕES E RESULTADOS

Este capítulo sintetiza de maneira sucinta os resultados obtidos ao longo deste projeto

de formatura. Quatro casos que possuem solução analítica, seis casos de uma mesma malha

de aterramento padrão e um caso de uma malha real de subestação foram ensaiados

utilizando-se a ferramenta de elementos finitos Ground 3D.

Os casos que possuem solução analítica foram simulados para que se pudesse validar o

programa. É vital que se obtenha resultados próximos aos analíticos nesta etapa para que o

modelo seja validado e para que possa haver o mínimo de confiança para simular casos

mais complexos que não possuem solução analítica, tomando o resultado do MEF como

próximo do esperado.

Uma malha padrão de 16 x 16 m foi ensada para seis perfis de solo diferentes. Os

resultados numéricos obtidos foram comparados com os previstos pela norma IEEE Std 80,

já que não há solução analítica para este caso. O solo foi estratificado em duas camadas e

cada caso considerava resistividades diferentes, para que se pudesse averiguar a

sensibilidade do resultado à estratificação.

O último caso analisado foi o da malha de aterramento de uma subestação de uma usina

hidroelétrica. Esta malha foi escolhida por ser representativa dos sistemas aterramentos de

subestações existentes no sistema elétrico brasileiro. Sua geometria e perfil do solo vieram

de plantas e memoriais de cálculos obtidos pelos autores junto ao lugar que trabalham.

Os onze casos simulados neste trabalho foram exaustivamente testados até que se

chegasse aos resultados aqui apresentados. É muito comum, durante uma campanha de

simulação, que os resultados sejam descartados por conta de discretizações não muito

adequadas, não convergência do método, travamento do computador por falta de memória

disponível, solução espúria sem significado físico, etc. Em razão disso, tomou-se bastante

tempo realizar estas simulações. Os resultados aqui apresentados foram os melhores que se

puderam obter com os recursos computacionais disponíveis.

45

5.1 COMPARAÇÃO COM MÉTODOS ANALÍTICOS

Os quatro casos de aterramento que possuem solução analítica, citados no capítulo 3,

foram simulados no GROUND 3D com o propósito de se validar o programa. Para todos os

casos, supôs-se que havia disponível 32 m de condutor de 4 mm de raio, a resistividade

aparente do solo considerada foi de 450 Ω.m e que os sistemas horizontais seriam

enterrados a 0,8 m de profundidade.

O primeiro caso analisado foi o de uma haste delgada de 32 m de comprimento,

enterrada verticalmente em um solo homogêneo, com uma corrente de defeito de 1 kA. De

acordo com a equação 3.1 esse aterramento deveria oferecer uma resistência de 20,98 Ω.

Ω=

−

⋅=

−

= 98,201004,0

324

322

4501

4

2n

r

Ln

LR a

h llππ

ρ

Equação 5-1

A haste foi simulada no Ground 3D, com uma malha de elementos finitos de 140 mil

nós (24 Mb) que levou cerca de cinco minutos para convergir. Neste caso, para diminuir o

tamanho da malha e, conseqüentemente, o tempo de convergência, tirou-se partido da

simetria do problema e simulou-se apenas um quadrante do domínio.

O MEF forneceu uma resistência de 21,84 Ω, muito próxima do valor teórico, com

desvio de apenas 4,1%. O mesmo ocorreu com a elevação de potencial. Como a resistência

é calculada a partir desta última, espera-se o mesmo desvio. A comparação é resumida na

Tabela 5-1.

Tabela 5-1 comparação de métodos para uma haste de 32 m de comprimento, ρa = 450 Ω.m e r = 4 mm

Analítico Numérico Diferença

Resistência 20,98 Ω 21,84 Ω 4,1%

Elevação de Potencial 20,98 kV 21,84 kV 4,1%

46

O caso da haste delgada também possui solução analítica para a distribuição de

potencial ao longo do solo a partir da distância da haste. A figura X mostra essa

distribuição e o gráfico da Ilustração 5-1 compara o resultado analítico e o numérico obtido

com o Ground 3D. O Método dos Elementos Finitos se mostrou muito capaz de calcular a

distribuição de potencial com precisão, mesmo próximo do condutor, aonde o gradiente

(campo elétrico) é muito intenso.

0

5000

10000

15000

20000

25000

0 5 10 15 20 25

Poten

cial Elétrico [V

]

Distância (m)

Comparação de resultados para uma

haste de 32 m de comprimento

Analítico

Ground 3D

Ilustração 5-1 – distribuição do potencial elétrico para o caso da haste delgada de 32 m, ρ = 450 Ω.m e r = 4 mm.

Em uma próxima etapa comparou-se o MEF com três casos mais elaborados que

possuem solução analítica. Para estes casos compararam-se apenas os valores das

resistências, por ser um parâmetro representativo da precisão da simulação.

Como já esclarecido, nestes casos pressupôs-se um comprimento total de condutores de

32 m, com raio de 4 mm, em solo de resistividade 450 Ω.m e enterrados a uma

profundidade de 0,8 m. O interessante em se fixar o comprimento de cabo é permitir uma

comparação entre os sistemas de aterramento.

47

Aplicando-se as equações 3.3, 3.4 e 3.5 para as geometrias já citadas, encontram-se as

resistências de terra das três configurações estudadas.

Para o condutor disposto horizontalmente:

Ω=

−+−

⋅= 00,24

32

8,022

8,0004,0

32

322

450 2

LlnRπ

Equação 5-2

Para o condutor disposto em L:

Ω=

++−

⋅= 78,24

16

8,04292,02373,0

8,0004,0

16

164

450 2

LlnRπ

Equação 5-3

Para o condutor disposto em cruz:

Ω=

+−+

⋅= 20,28

16

8,0284,4526,1

8,0004,0

16

164

450 2

LlnRπ

Equação 5-4

A Tabela 5-2 resume os resultados obtidos analiticamente e numericamente e estabelece

a comparação e a Ilustração 5-2 demonstra a distribuição de potencial ao longo do plano do

solo. Para os quatro casos estudados, o erro não ultrapassa 7,4%. Isso sugere que o MEF

sempre erra para cima, em razão das reflexões numéricas. Os valores dos desvios são

aceitáveis e os autores deste trabalho tomaram o programa GROUND-3D como validado,

ainda mais se considerando o ótimo desempenho da ferramenta ao calcular a distribuição de

potencial, mesmo nas proximidades dos condutores, onde a magnitude dos campos elétricos

é muito intensa.

Esta etapa foi muito importante porque os resultados aqui obtidos deram segurança aos

autores para que prosseguissem com simulações de geometrias mais complexas. Além

disso, os parâmetros PML, a discretização e tamanhos de domínio sugeridos pela referência

[3] se mostraram muito convenientes levando a resultados bem coerentes com tempos de

simulação relativamente baixos.

48

Tabela 5-2 – comparação de métodos para quatro casos analíticos

Configuração Tamanho do arquivo [Mb]

Tempo de simulação [min]

R [Ω] numérico

R [Ω] analítico

Diferença

(a) Haste vertical 24 5 21,84 20,98 4,1%

(b) Condutor enterrado 140 20 25,78 24,00 7,4%

(c) Condutores em L 195 20 26,52 24,78 7,0%

(d) Condutores em cruz 141 35 29,14 28,20 3,3%

(a) (b)

(c) (d)

Ilustração 5-2 – distribuição de potencial elétrico para quarto casos que possuem solução analítica

49

5.2 COMPARAÇÃO COM O MÉTODO DA NORMA IEEE STD 80

Uma malha de 16 x 16 m reticulada em 16 quadrados de 4 m de lado foi utilizada na

comparação do MEF com a norma IEEE Std 80 e sua geometria é dada pela Ilustração 5-3.

A malha não possui hastes em seus nós, é constituída por 160 m de condutor de 120 mm2 e

se situa a uma profundidade de 0,6 m.

Seus parâmetros para a fórmula de Sverak são: Lt = 160 m, Am = 256 m2, hm = 0,6.

16 m

16 m

Ilustração 5-3 – geometria da malha utilizada na comparação do MEF com a norma IEEE Std 80

Esta malha foi simulada para seis perfis de solo, estratificados em duas camadas de

resistividades distintas. A fórmula de Sverak é escrita em termos da resistividade aparante

na profundidade da malha e, por conta disso, ela foi determinada para uma profundidade p

a partir das duas camadas de cada perfil de solo através da fórmula de Schlumberger [2]:

( )

++=

21

1

21

21

ph

Ka ρρ Equação 5-5

Sendo:

12

12

ρρρρ

+−=K Equação 5-6

50

Onde:

K é o coeficiente de reflexão do solo

h1 é a profundidade da primeira cada [m]

p é a profundidade para a resistividade aparente [m]

ρ1 é a resisitividade da primeira camada [Ω.m]

ρ2 é a resisitividade da segunda camada [Ω.m]

Através da fórmula de Sverak, pode-se escrever a resistência da malha 16 x 16 m em

função da resistividade aparente na profundidade de 0,6 m:

aamR ρρ 032194,0256/206,01

11

25620

1

160

1 =

++

⋅+= Equação 5-7

As malhas 16 x 16 m nos seis perfis de solo foram simuladas numericamente e

calculadas pelo método proposto na norma IEEE Std 80. Os resultados são sumarizados na

Tabela 5-3. A malha de elemento finitos possui 400 mil nós (cerca de 48 Mb) e os tempos

de simulação foram de cerca de 10 minutos para cada caso.

Tabela 5-3 – comparação do MEF com a norma IEEE Std 80 para a malha 16 x 16 m

caso ρ1

[Ω.m]

h1

[m]

ρ2

[Ω.m] K

ρa

[Ω.m]

R [Ω] IEEE

Std 80

R [Ω]

numérico diferença

(a) 200 4,0 800 0,60 218 7,0 11,0 57%

(b) 3.000 4,0 100 -0,94 2.580 83,1 41,8 -50%

(c) 3.000 4,0 1.200 -0,43 2.808 90,4 61,7 -32%

(d) 3.000 4,0 10.000 0,54 3.232 104,4 151,9 46%

(e) 3.000 1,2 100 -0,94 1.639 52,8 20,9 -60%

(f) 3.000 0,6 100 -0,94 490 15,8 9,2 -41%

Há uma grande discrepância entre o resultado numérico e o calculado pela norma,

porém os valores se situam na mesma ordem de grandeza. Isto já era esperado pelo fato do

modelo da norma ser semi-empírico, isto é, não é exato, pois há muitas hipóteses

simplificadoras. Para estas simulações, a distribuição de potencial elétrico está mostrada na

Ilustração 5-4.

51

(a) (d)

(b) (e)

(c) (f)

Ilustração 5-4– distribuição de potencial para os seis perfis de solo para a malha 16 x 16 m

52

5.3 COMPARAÇÃO COM UMA MALHA REAL

Com o fim de se ensaiar o método dos elementos finitos como ferramenta de projeto de

malhas maiores e mais complexas, como as encontradas nos sistemas de potência, os

autores deste trabalho simularam a malha da subestação do pátio de 230 kV de uma

importante usina hidroelétrica brasileira, cujos dados foram obtidos mediante autorização

da empresa em que trabalham os autores.

A malha dessa usina é irregular, tem 81 m de comprimento por 49,5 m de largura, sendo

os condutores espaçados de maneira não-simétrica. Sua geometria é dada pela Ilustração

5-5, as dimensões estão todas em metros.

81 m

49,5 m

11,5

1

1,5

12,

0

11,

5

3,0

16,0 17,0 15,0 16,0 17,0

Ilustração 5-5 - geometria da malha de aterramento simulada neste trabalho

O memorial de cálculo desta malha se utiliza da norma IEEE Std 80 para estimar o

valor de sua resistência de terra. Mesmo não sendo regular e simétrica, a malha satisfaz as

hipóteses e premissas da fórmula de Sverak viabilizando a utilização do método. Os

parâmetros da fórmula para esta geometria são: Lt = 783 m, Am = 4.009,5 m2, hm = 0,3. A

resistividade aparante do solo, para a profundidade da malha é de ρa = 2.113 Ω.m.

53

Logo, para esta malha, a resistência da tomada de terra pelo IEEE Std 80 vale:

Ω=

++

⋅+= 46,17

5,4009/203,01

11

5,400920

1

783

12113mR Equação 5-8

Quanto à simulação numérica, inicialmente surgiram muitas dificuldades em virtude do

tamanho da malha de elementos finitos. A malha com a discretização típica chega a ter

quatro milhões de nós e ocupar 630 Mb de memória. Em razão da impossibilidade de se

simular um caso tão grande, houve a necessidade de se reduzir a discretização do domínio

com o aumento do tamanho dos elementos, levando a 1,6 milhões de nós. Entretanto, o

baixo valor do potencial elétrico calculado na fronteira que representa o ponto remoto

garantiu que o resultado final era coerente, mesmo com uma malha menos refinada. Para

esta simulação utilizou-se um computador dual-core de 2,8 GHz e com 2 GB de memória

RAM. A malha final resultou em um arquivo de nós e elementos de 245 Mb e levou 1 hora

e 55 minutos para atingir a convergência.

A resistência de aterramento numérica resultou em 8,45 Ω frente a 17,46 Ω previsto

pelo IEEE Std 80. Novamente, o fato é explicado pelas hipóteses simplificadoras da norma.

Os resultados foram comparados na Tabela 5-4.

Tabela 5-4 – comparação de métodos para a malha de uma subestação típica, ρa = 450 Ω.m e r = 4 mm. número de elementos

tamanho da malha

tempo de simulação

numérico IEEE Std 80 diferença

1,6 milhões 245 Mb 115 min 8,45 Ω 17,46 Ω 51,6%

Quanto à distribuição do potencial elétrico, através da Ilustração 5-6 percebe-se que

a malha cumpre bem o seu papel de equalização de potencial, ao menos dentro da

subestação, minimizando os riscos a um operador transeunte quando do surgimento de uma

falta a terra. Porém, nas extremidades da malha há um grande gradiente de potencial, que

corresponde a campos elétricos intensos. O pior caso é ilustrado na Ilustração 5-7, onde se

traça o perfil de potencial elétrico ao longo da diagonal de uma extremidade da malha para

o caso de uma corrente de defeito de 25 kA. Nota-se que há pontos ao longo da diagonal

aonde o gradiente chega a quase 4 kV/m.

54

Ilustração 5-6 – distribuição de potencial elétrico na malha

180

185

190

195

200

205

0 5 10 15 20 25

Potencial Elétrico [kV]

distância percorrida na diagonal sobre o solo [m]

Distribuição de potencial entre dois condutores para uma corrente de defeito de 25 kA

Ilustração 5-7 – distribuição de potencial ao longo de uma diagonal da malha

55

5.4 RESUMO DAS COMPARAÇÕES

Neste trabalho, foram simulados quatro casos de aterramento que possuem solução

analítica, seis malhas de terra regulares para comparação com a norma IEEE Std 80 e um

caso real, a malha da subestação de uma usina hidroelétrica.

A utilização do MEF nos casos que possuem solução analítica forneceu os subsídios

para a validação do programa Ground 3D, já que os resultados obtidos numericamente

foram muito próximos daqueles analíticos, tanto para a distribuição de potencial elétrico

quanto para a resistência de terra.

Em uma próxima etapa, utilizou-se de uma malha regular, de 16 x 16 m, para avaliar o

desempenho do MEF na simulação de malhas de terra e compararam-se os resultados

obtidos com aqueles calculados através da norma IEEE Std 80. Como o método da norma é

semi-empírico, oriundo de muitas simplificações e hipóteses como a resistividade aparente,

tomam-se esses valores com cautela e já se espera que eles sejam apenas uma previsão. O

desvio do resultado numérico do previsto pela norma é de cerca de 50%, em média.

Admitindo-se que resultado numérico seja provavelmente mais próximo do real do que

o método aproximado da norma, aparentemente o método da norma erra para cima quando

a resistividade da primeira camada é maior do que a da segunda (K negativos) e para baixo

quando a resistividade da primeira camada é menor do que a da segunda (K positivos). Isso

sugere que a fórmula de Sverak subestima a influência da segunda camada, porém esse fato

não pode ser afirmado sem que mais investigações sejam feitas.

Por fim, simulou-se a malha da subestação de uma usina hidroelétrica com o intuito de

se testar a ferramenta de elementos finitos (Ground 3D) em uma malha significativamente

maior e mais complexa e, por conta disso, surgiram inúmeras dificuldades. O primeiro

empecilho foi o tamanho da malha de elementos finitos, que consumia praticamente toda a

memória do computador. A solução foi simular em um computador pessoal de alto

desempenho, porém, mesmo assim, os tamanhos de malhas ficavam intratáveis quando se

utilizava a discretização proposta por [3]. Neste caso, houve a necessidade de se engrossar

um pouco a malha, perdendo ligeiramente em precisão, porém os resultados obtidos foram

aceitáveis e coerentes.

56

Outro problema que surgiu foi a dificuldade na manipulação dos arquivos que serviam

como dados de entrada. Os autores tiveram que desenvolver programas computacionais em

linguagem C para gerar arquivos de condição de contorno, arquivos de parâmetros PML,

etc. Mais tarde, parte desses problemas foi resolvida ao se utilizar uma versão mais nova do

Ground 3D, que possui uma melhor interface com o GMSH. Também foi elaborada uma

planilha em Excel para que se pudesse ler o valor do potencial elétrico nos nós da malha e

gerar gráficos de distribuição de potencial ao longo de caminhos traçados no solo.

Superadas as dificuldades, foi possível simular a malha da usina em um computador

pessoal de alto desempenho, com tempo convergência de 1 hora e 55 minutos. Os

resultados obtidos demonstraram ser coerentes, tanto pela inspeção visual da distribuição de

potencial elétrico, quando pelo seu valor na fronteira, aonde foi bem amortecido. Além

disso, o valor da resistência de terra está em uma faixa condizente com o esperado. Julgou-

se, então, que a simulação foi válida e representativa.

57

6. CONCLUSÕES

Há mais de vinte anos, o Laboratório de Eletromagnetismo Aplicado Prof. Dr. Áurio

Gilberto Falcone (LMAG) tem sido pioneiro na aplicação do método dos elementos finitos

(MEF) em sistemas de aterramento. Este laboratório desenvolveu inúmeras ferramentas

computacionais para estudos de eletromagnetismo, em particular o GROUND 3D, que foi

utilizado neste trabalho. Porém, as malhas de aterramento estudadas até então se tratavam,

em sua maioria, de casos acadêmicos. Este trabalho teve um viés prático, voltado para as

malhas de aterramento.

Ainda que se tenha feito muitos progressos no desenvolvimento destas ferramentas e a

capacidade computacional crescido consideravelmente neste par de décadas, pouquíssimos

são os casos de malhas de aterramento de subestações que foram projetadas por métodos

numéricos. Ainda hoje frequentemente se utiliza a norma IEEE Std 80, que se baseia em

modelos semi-empíricos, para este tipo de projeto.

Contudo, com o advento das subestações compactas com equipamentos isolados a SF6 e

a tendência de se trazer as subestações para próximo dos consumidores, em meio a grandes

cidades, soluções não convencionais para a dissipação da corrente de falta têm sido

adotadas. Por essa razão espera-se que o método numérico seja cada vez mais utilizado,

face à inexistência de uma metodologia de projeto adequada para estes casos que não a

experimental. Esses fatos foram a motivação necessária para que os autores prosseguissem

seu trabalho neste tema.

A comparação de resultados com os casos que possuem solução analítica provaram que

o MEF, em particular sob a forma do programa GROUND 3D, é uma ferramenta muito

capaz na previsão acurada da distribuição de potencial elétrico e na estimação da resistência

a terra. A confiança adquirida em razão dos excelentes resultados obtidos nesta etapa deu

subsídio para que análises mais complexas pudessem ser feitas.

58

Em um segundo momento, comparou-se a resistência prevista pela norma IEEE Std 80

com os resultados numéricos de uma malha simples, reticulada, para diversas configurações

de solo. Como não há solução analítica para uma malha de terra, confiou-se no resultado

numérico e admitiu-se que ele seria muito próximo da solução exata, pois a inspeção visual

da distribuição de campos se mostrou muito coerente, o método convergiu e o potencial no

ponto remoto era sempre muito pequeno. Logo se concluiu que a previsão da fórmula de

Sverak, que é o método da norma já citada, é um pouco distante do resultado exato e isso já

era esperado por conta da natureza semi-empírica do método. Porém, ainda assim a fórmula

de Sverak provê uma primeira estimativa da resistência de uma malha de aterramento, pelo

menos da sua ordem de grandeza.

Talvez esse resultado justifique a utilização quase que unânime da norma IEEE Std 80

no projeto de subestações. Soma-se a isso o fato da resistividade aparente do solo ser uma

variável ainda mais sujeita a erros, sendo uma função da umidade, da temperatura e do

tempo, e então o projeto de uma malha passa a ser praticamente heurístico.

Porém, este método semi-empírico não provê a distribuição de potencial ao longo do

solo, apenas algumas fórmulas semi-empíricas para estimar os potenciais máximos de toque

e passo presentes, que também podem não ser muito precisos. Por sua vez, um método

numérico, como o MEF, possibilita a visualização dessa distribuição e os potenciais

máximos podem ser estimados com muito mais precisão.

Por conta dessas vantagens intrínsecas aos métodos numéricos, os autores propuseram a

simulação de uma malha de aterramento real, representativa daquelas existentes em

qualquer grande sistema de potência. Dessa forma, foi possível analisar o desempenho do

MEF, em particular do GROUND 3D, como ferramenta de projeto de malhas maiores e

mais complexas.

No início houve muita dificuldade em gerar uma malha com uma discretização

adequada à simulação. O tamanho do domínio do estudo frente ao tamanho da malha de

terra e de seus condutores fez com que o número de elementos necessários aumentasse

significativamente, a ponto de exigir muita memória do computador e, por muitas vezes,

causando a instabilidade do sistema operacional. Mesmo utilizando-se um computador

pessoal de alto desempenho, foi necessário que a malha desse caso fosse engrossada um

pouco para que houvesse uma diminuição do número de elementos.

59

Boa parte dos problemas de falta de memória poderiam ter sido resolvidos se os autores

tivessem utilizado outros sistemas operacionais baseados em Linux ou Unix ao invés do

Windows XP. Estes outros sistemas possuem um gerenciamento de memória muito

superior, podendo alocar memória de maneira mais eficiente, além deles próprios ocuparem

pouca memória RAM. Algumas estações de trabalho de alto desempenho do LMAG já

possuem esses sistemas alternativos.

Surgiram muitas outras dificuldades, principalmente com os arquivos de dados, e os

autores foram obrigados a construir programas para automatizar a manipulação desses

arquivos. Superados os problemas, a simulação ocorreu de forma tranqüila, demorando

cerca de duas horas para atingir a convergência. Os resultados obtidos foram muito

satisfatórios e coerentes com o que era esperado. Conseguiu-se traçar a distribuição de

potencial ao longo da diagonal de um dos reticulados da malha, o mais próximo à borda,

aonde os campos elétricos eram mais intensos. Esse tipo de traçado pode ser utilizado para

estimar com maior precisão o potencial de passo da malha. Fica como sugestão de trabalho

futuro utilizar o MEF para comparar os potenciais de passo numéricos com aqueles obtidos

através da norma e de medições em campo, se for aplicável.

Ainda assim, o projeto de um sistema de aterramento através do MEF deve ser cercado

de uma série de cuidados. Principalmente, a malha deve ser muito bem preparada,

discretizada convenientemente para que não fique muito pesada e ainda assim permita uma

precisão adequada. A utilização da ferramenta per se não garante bons resultados. Cabe ao

usuário discernimento para interpretar os dados e julgar a qualidade do resultado.

O MEF se mostrou uma ferramenta muito capaz na resolução de problemas

eletromagnéticos, notadamente os de aterramento, e deve ser uma ferramenta muito

utilizada em projetos futuros de subestações, mais complexas, em meio a grandes cidades,

com aterramentos naturais próximos, como encanamentos, além da possibilidade de haver

potenciais transferidos de outras malhas próximas. Tudo isso é passível de ser analisado em

um programa como o GROUND 3D, bastando apenas pequenas modificações em seu

código fonte, caso seja necessário.

60

É importante lembrar que ao longo destes últimos vinte anos, muitos avanços foram

feitos sobre o programa original. Devido ao trabalho de muitas pessoas, surgiram versões

mais modernas, como a implantação do MEF de arestas, melhorias na biblioteca de solução

do sistema linear, a abordagem do elemento unidimensional com salto de potencial, a

utilização do PML para redução do erro de truncamento, a alimentação do sistema em

corrente ao invés de tensão, o estudo do problema de aterramento em regime harmônico

através da formulação A-V, a utilização do GMSH para geração da geometria, da malha de

elementos finitos e do pós-processamento, etc. Há muito campo para desenvolvimento

nesta área e num futuro próximo, como o avanço dos sistemas computacionais, as

possibilidades serão praticamente ilimitadas.

Ao longo deste último par de décadas, o amadurecimento dos computadores e do

programa de elementos finitos para aterramento implantado no LMAG viabilizou a

simulação de uma grande malha de subestação, o que não era possível naqueles primeiros

anos de desenvolvimento. Este trabalho atingiu o seu objetivo ao demonstrar que, no estado

atual, é possível o estudo de complexos sistemas de aterramento em regime estacionário

através do MEF a partir de modestos recursos computacionais.

61

7. BIBLIOGRAFIA

[1] G. Kindermann, J. M. Campagnolo. Aterramento Elétrico 4º edição, Sagra Luzatto, 1998, Porto Alegre.

[2] ANSI/IEEE Std 80 – 2000 – IEEE Guide for safety in AC substation grounding – Pub. by IEEE. New York, John Wiley, 2000.

[3] SILVA, V. C. Método dos Elementos Finitos aplicado à solução de problemas de aterramento elétrico. São Paulo, 2006. 89p. Tese (Livre Docência) – Escola Politécnica da Universidade de São Paulo.

[4] J. R. Cardoso, Ground-3D: Uma Contribuição à Análise dos Sistemas de Aterramento pelo Método dos Elementos Finitos. São Paulo, 1993. 71p. Tese (Livre Docência) – Escola Politécnica da Universidade de São Paulo.

[5] J. R. Cardoso, Problemas de campos eletromagnéticos estáticos e dinâmicos : uma abordagem pelo método dos elementos finitos, 1985. 159p. Tese (Doutorado) – Escola Politécnica da Universidade de São Paulo.

[6] J. R. Cardoso, Introdução ao Método dos Elementos Finitos para Engenheiros Eletricistas, Publicação Independente, 1ª Edição, 1995, São Paulo.

[7] Ch. Geuzaine, J. F. Remacle, Gmsh© A three-dimensional finite element mesh generator with built-in pre-and-post-processing facilities, Version 2.0.8, Agosto 2007, http://www.geuz.org/gmsh/

[8] J. Jin, The Finite Element Method in Electromagnetics, 2nd edition, 2002, John Wiley & Sons Inc. NY.

[9] J. P. Bérenger, “A perfectly matched layer for the absorption of electromagnetic waves”, J. Comp. Phys., Vol. 114, Oct. 1994, 185-200.

[10] P. P. Silvester, R. L. Ferrari, Finite Elements for Electrical Engineers, Cambridge University Press, 1996, 3ª edição.

[11] Matthew N. O. Sadiku, Numerical Techniques in Electromagnetics, Boca Raton, CRC Press, 1992