Campo e Potencial Elétrico numa Região do Espaço Confinada por uma
Fronteira Fractal.
Thiago Albuquerque de Assis
Campo e Potencial Elétrico numa
Região do Espaço Confinada por
uma Fronteira Fractal.
Thiago Albuquerque de Assis
Orientador: Prof. Dr. Caio Mário Castro de Castilho
Co-Orientadores: Prof. Dr. Fernando de Brito Mota
Prof. Dr. José Garcia Vivas Miranda
Dissertação apresentada ao Instituto de Fı́sica da Universidade Federal da Bahia
como requisito parcial para a obtenção do grau de Mestre em Fı́sica.
Março de 2007
Aos meus queridos pais Denise e Roberto e ao meu irmão Raphael.
Agradecimentos
A Deus pela sua presença constante na minha vida, me proporcionando esta
magnı́fica oportunidade de poder estudar ciências.
À Pós-Graduação do IF pela preocupação notável na minha formação, bem
como na realização deste estudo.
Ao Professor Sérgio Cavalcante Esperidião (In Memorian), por ter me incentivado nos estudos e objetivos.
Ao Colega Hugo de Oliveira Dias Filho (In Memorian), por ter me familiarizado
com o problema a ser solucionado no inı́cio da minha carreira cientı́fica.
Aos Professores David Vianna e Caio Castilho pela formação adquirida durante
a realização deste trabalho.
Aos Professores Caio Castilho, Fernando de Brito, José Garcia, Roberto Andrade, Hélio Silva Campos e Anderson Barbosa que viram neste trabalho uma proposta
significativa de tratar fenômenos naturais e pelo relevante auxı́lio na produção desta
dissertação.
Ao Grupo de Fı́sica de Superfı́cies e Materiais por contribuições significativas
na minha formação acadêmica.
À minha namorada Viviane, pela paciência e compreensão em situações que
precisei abdicar dos momentos de lazer.
Ao Professor Paulo Bahiense por me ensinar a ensinar.
Ao Amigo Frei Paulo Avelino por buscar sempre me incentivar nos momentos
em que precisei de auxı́lio, durante este trabalho.
A todos os meus amigos do mestrado que souberam compartilhar comigo,
momentos de alegria e tristeza.
Resumo
Neste trabalho considera-se o comportamento do campo elétrico numa região
delimitada por um perfil ou por uma superfı́cie com geometrias fractais, e uma linha
ou um plano distante. Entre as fronteiras é mantida uma diferença de potencial
constante, assumindo cada uma das fronteiras como possuindo propriedades condutoras. A solução da equação de Laplace foi numericamente obtida pelo método de
Liebmann e os expoentes de rugosidade, bem como as dimensões fractais das linhas
equipotenciais, foram calculados através dos métodos RMS e Semivariograma. Os
resultados indicam, para os cálculos desenvolvidos no caso 1D+1, que para uma
determinada linha equipotencial situada entre as fronteiras, a dimensão fractal cresce
com o tamanho do comprimento do perfil, sugerindo que todas as linhas equipotenciais
devem ser caracterizadas pelo mesmo valor de dimensão fractal, para sistemas com
comprimento infinito. Calcularam-se também as dimensões fractais das superfı́cies
equipotenciais para condutores base com diferentes geometrias. Estas últimas apresentam uma rugosidade que decresce à medida que a distância média da correspondente
superfı́cie equipotencial ao condutor base cresce. Uma vez que as micro-irregularidades
do condutor base se refletem no comportamento do campo elétrico local, elas afetam as
trajetórias de partı́culas carregadas não interagentes sujeitas a este campo. Utilizando
os métodos de Dinâmica Molecular Clássica, discute-se as “informações”, contidas nas
trajetórias destas partı́culas e nas suas grandezas dinâmicas, a respeito da morfologia
da superfı́cie irregular.
Parte do conteúdo do Capı́tulo 3 corresponde a material publicado em Journal
of Physics: Condensed Matter, 18, 3393 (2006). O conteúdo do Capı́tulo 4 deverá ser
submetido a publicação.
v
Abstract
In this work it was considered the behavior of the electric field in a region
bounded by a linear profile, or by a surface, with a fractal geometry and a distant
straight line or a plane. Between the boundaries it is kept a constant potential
bias, assuming both boundaries as being conducting ones. The solution for Laplace’s
equation was numerically obtained by Liebmann’s method while the rough exponents
and fractal dimensions were calculated using the RMS and Semivariogram methods.
The results indicate, for the performed calculations in 1D+1, that for a fixed equipotential line between the boundaries, its fractal dimension increases with the profile size,
suggesting that all the equipotential lines should be characterized by the same fractal
dimension value, in the case of systems with an infinite size. It was also calculated
the fractal dimensions for equipotential surfaces in the case of basal conductors with
different geometries. These last present a roughness that decreases with the increase
of the average distance from the equipotential to the basal conductor. Since the
micro-irregularities of the basal conductor are reflected in the behavior of the local
electric field, they affect the trajectories of non-interagent particles under the action of
the field. Using the methodology of Classical Molecular Dynamics, it is explored the
“information”, contained within the particles trajectories and dynamical quantities,
concerning the irregular surface morphology.
Part of the content of Chapter 3 corresponds to a material published in Journal
of Physics: Condensed Matter, 18, 3393 (2006). The content of Chapter 4 will be
submitted to publication.
vi
Conteúdo
1 Introdução
1
2 Alguns Aspectos a Respeito da Geometria Fractal
5
2.1
Um Breve Histórico . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2
Dimensão Fractal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3
Fractais Auto-Similares e Fractais Auto-Afins. . . . . . . . . . . . . . . 14
2.4
Métodos para Determinação da Dimensão Fractal . . . . . . . . . . . . 17
2.5
5
2.4.1
Homotetia e Dimensão Fractal . . . . . . . . . . . . . . . . . . . 17
2.4.2
Método de Contagem de Caixas. . . . . . . . . . . . . . . . . . . 19
2.4.3
Método RMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.4.4
Método Semivariograma . . . . . . . . . . . . . . . . . . . . . . 25
2.4.5
Análise de Fourier
. . . . . . . . . . . . . . . . . . . . . . . . . 28
Simulação Computacional de Superfı́cies Fractais . . . . . . . . . . . . 29
2.5.1
Técnica FBM . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.5.2
Deposição Balı́stica . . . . . . . . . . . . . . . . . . . . . . . . . 33
3 Campo Elétrico Gerado por Fronteiras Auto-Afins
39
3.1
Equações Diferenciais Parciais . . . . . . . . . . . . . . . . . . . . . . . 40
3.2
Solução Analı́tica da Equação de Laplace . . . . . . . . . . . . . . . . . 41
3.3
Solução Numérica em Duas Dimensões . . . . . . . . . . . . . . . . . . 44
3.4
Discussão e Resultados (1D+1) . . . . . . . . . . . . . . . . . . . . . . 45
3.5
Discussão e Resultados (2D+1) . . . . . . . . . . . . . . . . . . . . . . 55
vii
4 Dinâmica de Partı́culas Carregadas em um Campo Elétrico Irregular 70
4.1
Equações do Movimento . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.1.1
4.2
Método Preditor-Corretor . . . . . . . . . . . . . . . . . . . . . 72
Discussão e Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.2.1
Estudo do Campo Elétrico Médio . . . . . . . . . . . . . . . . . 77
4.2.2
Movimento de Partı́culas não Interagentes em um Campo Elétrico
Irregular . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5 Conclusões e Perspectivas
99
A Caracterı́sticas e Propriedades de Alguns Fractais Auto-Similares 106
A.1 Conjunto de Cantor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
A.2 Ilha de Koch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
A.3 Triângulo de Sierpinski . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
A.4 Esponja de Menger . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
B Equação KPZ
116
C Método de Liebmann
119
D Algoritmo de Verlet
125
viii
Capı́tulo 1
Introdução
Atualmente ainda existem muitos problemas relacionados à Eletrostática Clássica cuja solução não é trivial. Entre estes está incluı́da a solução analı́tica das
equações de Poisson e de Laplace em regiões onde a fronteira não é regular. Um outro
exemplo, este de utilidade prática, diz respeito ao cálculo da emissão por campo,
seja no caso de elétrons que, por tunelamento, são emitidos a partir de extremidades
condutoras ou semicondutoras, como também na emissão de ı́ons produzidos pela
transferência de elétrons a partir de átomos ou moléculas. Uma questão importante
nestes fenômenos diz respeito à interação do campo eletrostático com a superfı́cie
do emissor, cuja geometria muitas vezes é irregular. Nestes casos esta interação não
possui uma solução analı́tica e, portanto, métodos numéricos constituem a alternativa
a ser utilizada. Ignorar a exata morfologia da estrutura emissora, substituindo-a por
uma mais simples no que diz respeito ao tratamento matemático, é uma prática
comum. Isto muitas vezes permite a determinação da magnitude do campo eletrostático produzido em uma superfı́cie condutora, e daı́ são estimados parâmetros
geométricos para os dispositivos emissores, estes às vezes possuindo uma geometria
idealizada, onde são eliminadas rugosidades. Tais relações analı́ticas acabam por
serem capazes de tratar geometrias como elipsóides, hiperbolóides superimpostos a
planos ou outras formas de suporte; enfim, um conjunto vasto de emissores cujas
geometrias apresentam protuberâncias que representam as irregularidades caracterı́sticas da morfologia de emissores em escala nanométrica. A possibilidade de tratar
1
superfı́cies irregulares e os campos na vizinhança destas resulta assim num passo
adiante nas aproximações utilizadas.
Entender o campo eletrostático e sua interação com a matéria é uma questão
de fundamental importância em vários dispositivos.
Isto pode ser exemplificado
fazendo-se uma referência à microeletrônica a vácuo, onde o cuidado com o cálculo da
distribuição do campo elétrico próximo à superfı́cie emissora é essencial, possibilitando
o cálculo da área de emissão e a corrente total de emissão se a função trabalho é
conhecida. A própria determinação experimental da função trabalho requer uma
adequada consideração a respeito do campo existente na região vizinha à superfı́cie.
Como é possı́vel constatar [1] que a emissão ocorre, essencialmente, através de irregularidades em pequenas escalas, resultado do aumento local da intensidade do campo
elétrico, torna-se importante a determinação de parâmetros geométricos para os emissores. Isto decorre, em última instância, do fato que, experimentalmente, tem-se
controle sobre a tensão e a intensidade de corrente, não sobre o campo ou a densidade
local de corrente, estes últimos os fatores usualmente calculados a partir de uma
abordagem analı́tica.
Na análise da morfologia de uma superfı́cie torna-se essencial o conceito de
escala, usada pela Mecânica Estatı́stica moderna para demonstrar os chamados comportamentos universais de escala, ou seja, para mostrar que sistemas, aparentemente
diferentes, apresentam um comportamento de escala comum. Existem, portanto,
certas “leis de escala” básicas e independentes de muitos detalhes desses sistemas.
Logo, a caracterização de sistemas através de expoentes globais leva à definição de
classes de universalidade. Dois sistemas pertencem à mesma classe universal se podem
ser descritos pelos mesmos expoentes de escala [2]. Em outras palavras, este tipo de
descrição é uma tentativa de evidenciar novas ”simetrias” embutidas em sistemas
aparentemente distintos e sem qualquer lei de formação.
De acordo com dados experimentais, as superfı́cies de alguns materiais apresentam auto-similaridades em certas escalas e, portanto, podem ser descritas por modelos
fractais [3] . A geometria fractal, bem como os fenômenos caóticos, têm sido, nos
últimos anos, alvo de contı́nuas investigações em todo o mundo. As técnicas associadas
2
à geometria fractal têm-se revelado uma ferramenta extremamente útil a muitas
ciências. Estas, que de uma forma ou de outra tratam de interfaces num universo
despovoado de formas geométricas perfeitas, estudam fenômenos onde proliferam
superfı́cies irregulares. Portanto, a geometria fractal apresenta-se como um meio
de tratar alguns fenômenos até agora considerados como erráticos, imprevisı́veis e
aleatórios. Para citar alguns exemplos, existem modelos fractais para tratamento
da deposição de filmes finos por MBE (Molecular Beam Epitaxy), por Vaporização
Catódica (Sputtering) [4] e outras técnicas. O interesse por fractais nessa área foi
impulsionado ainda pela possibilidade de se analisar superfı́cies utilizando microscópios
de alta resolução, como AFM (Atomic Force Microscope) e STM (Scanning Tunnelig
Microscope) [5], capazes de identificar formações em escala atômica.
No caso de superfı́cies, um conceito relevante é a medida de sua rugosidade.
Neste trabalho foram estudadas as proriedades do campo elétrico gerado por perfis e
superfı́cies com geometria fractal. Para o caso de um perfil fractal, foi verificado como
o tamanho do mesmo afeta as irregularidades das linhas equipotenciais, considerandose uma diferença de potencial constante entre o mesmo e uma linha reta distante. Para
o caso de superfı́cies, foram estudadas as superfı́cies geradas por deposição balı́stica e
as superfı́cies formadas iterativamente a partir de uma metodologia FBM (Fractional
Brownian Motion), através de um processo randômico auto-afim. Tratando estas
superfı́cies como condutoras e, atribuindo-se uma diferença de potencial elétrico entre
cada uma destas superfı́cies e um plano distante, foram determinadas também, as
propriedades do campo elétrico a partir da solução numérica da equação de Laplace
na região delimitada. Constatou-se, a partir da determinação das dimensões fractais
das superfı́cies equipotenciais que estas possuem comportamento de invariância de
escala e apresentam uma rugosidade que decresce à medida que cresce a distância
média da correspondente superfı́cie equipotencial ao condutor base. Uma vez que as
irregularidades do condutor base se refletem no comportamento do campo elétrico
local, verificou-se que as trajetórias de partı́culas carregadas, determinadas pelo método de Dinâmica Molecular Clássico, trazem informações relevantes a respeito da
morfologia de superfı́cies, sendo portanto de interesse prático para, por exemplo, a
3
determinação do comportamento de eletrodos [6].
Esta dissertação contém mais 4 capı́tulos. No Capı́tulo 2 é apresentada uma
discussão a respeito das principais definições associadas à geometria fractal. Métodos
para a determinação da dimensão fractal de perfis, bem como de superfı́cies também
são apresentados. Por fim, é feita uma abordagem a respeito da simulação computacional de superfı́cies fractais (Deposição Balı́stica (DB) e metodologia Fractional
Brownian Motion (FBM)). No Capı́tulo 3 é feita preliminarmente uma discussão a
respeito das equações diferenciais parciais, bem como das suas soluções, analı́tica e
numérica, em duas dimensões. Particular menção é feita para a solução da equação
de Laplace. São apresentados os resultados obtidos para o estudo de perfis bem como
de superfı́cies equipotenciais, além dos domı́nios de validade dos métodos utilizados
para o cálculo dos expoentes de rugosidade (RMS e Semivariogramna). No Capı́tulo
4 apresenta-se uma breve discussão a respeito do método de Dinâmica Molecular,
utilizando-se o método Preditor-Corretor. São apresentados os resultados obtidos na
investigação do comportamento médio da intensidade do campo elétrico gerado por
superfı́cies fractais com tamanhos 100x100. Por fim, são apresentados os resultados
obtidos via dinâmica de partı́culas carregadas em um campo gerado por superfı́cies
fractais. No Capı́tulo 5 são apresentadas as conclusões, discutindo-se perspectivas
para estudos futuros.
Este trabalho também contém quatro apêndices. No Apêndice A são apresentadas algumas caracterı́sticas e propriedades de alguns fractais auto-similares
(Conjunto de Cantor, Ilha de Koch, Triângulo de Sierpinski e Esponja de Menger).
No Apêndice B é apresentada a equação KP Z, que se apresenta como classe de
universalidade descrevendo matematicamente o modelo de deposição balı́stica em
largas escalas. No Apêndice C é discutido o método de Liebmann para a solução
numérica da equação de Laplace em 1D+1 e sua extensão para 2D+1 dimensões.
Finalmente, no Apêndice D, é apresentado o algoritmo de Verlet complementando
a discussão iniciada no Capı́tulo 4 com o algoritmo Preditor-Corretor.
4
Capı́tulo 2
Alguns Aspectos a Respeito da
Geometria Fractal
No presente Capı́tulo, será apresentada uma breve discussão a respeito de
alguns conceitos relativos ao que se denomina de Geometria Fractal. Além de um
brevı́ssimo histórico a respeito do tema, a definição formal de dimensão fractal e
uma descrição de métodos utilizados, ao longo deste trabalho, para seu cálculo, será
exposta. Por fim, dois métodos para a simulação computacional de superfı́cies com
geometria fractal serão descritos.
2.1
Um Breve Histórico
O emprego do termo fractal pode ser temporalmente localizado no ano de
1975, quando Benoit Mandelbrot pela primeira vez fez uso deste termo [7]. Há porém
uma história anterior, quando uma série de considerações e estudos, que sem os seus
protagonistas o soubessem, abriram caminho para o desenvolvimento do tema [7].
Entre a segunda metade do século XIX e a primeira metade do século XX,
foram propostos vários objetos matemáticos, com caracterı́sticas especiais, que foram,
durante muito tempo, considerados “monstros matemáticos”, uma vez que desafiavam
as noções comuns de infinito e para os quais não havia uma explicação objetiva.
Cantor (1845-1918), que se destacou por idéias altamente inovadoras sobre o conceito
5
de infinito, propôs a construção de um objeto que resultou chamar-se de poeira de
Cantor [8]. A proposta, formulada em 1895, era a de considerar um segmento linear,
do qual se removeria o seu terço médio. A partir daı́, o terço médio de cada um
dos dois segmentos restantes também seria removido, e assim sucessivamente com
cada porção resultante da remoção de terços médios (ver Figura 2.1) . A repetição
infinita do processo gera uma “poeira” que, sendo constituı́da de um número infinito
de pontos, ainda assim possui um comprimento total igual a zero. Sendo o Conjunto
de Cantor denotado por C e assumindo-o como contido no intervalo [0,1], define-se o
mesmo como sendo:
C=
∞
\
In
(2.1)
n=0
onde:
I0 ⊃ I1 ⊃ I2 ⊃ ... ⊃ In
(2.2)
Figura 2.1: Os três primeiros nı́veis para a construção do Conjunto de Cantor.
De modo similar, foi apresentada em 1904 a Curva de von Koch [9], que
corresponde a uma linha, de comprimento infinito, delimitando uma área finita. A
construção se dá a partir de um segmento de reta que em seguida é dividido em três
segmentos iguais. Depois disto, substitui-se o terço médio por um triângulo equilátero
tirando-lhe a sua base. O processo iterativo consiste em aplicar a mesma regra a
6
cada um dos segmentos de reta que resultam da iteração imediatamente anterior.
Se considerar-se cada passo, nota-se que, de um nı́vel para a seguinte, substitui-se
três segmentos por quatro de igual comprimento, ou seja, o comprimento total é
multiplicado por 4/3 em nı́veis sucessivos. O limite de uma sucessão geométrica de
razão 4/3 é infinito, o que significa que a figura final, i.e., aquela para a qual tende
a sucessão descrita, terá um comprimento infinito. Este limite foi designado por
Mandelbrot como “infinito interno”. Portanto, no n-ésimo nı́vel, o comprimento da
curva de Koch ser dado por:
Sn = Sn−1 +
Sn
4
= ( )n
3
3
(2.3)
No limite, para um número de infinitos nı́veis, tem-se o seguinte resultado:
lim Sn = ∞
n→∞
(2.4)
Esta tendência, de uma propriedade geométrica tender a infinito, é uma caracterı́stica, tı́pica dos fractais. Os seis primeiros nı́veis da Curva de Koch estão
mostrados na Figura 2.2. Observando o último nı́vel representado, tem-se a impressão
visual de que a curva parece ter uma certa “espessura”, devido às constantes mudanças
de direção no processo construtivo. Isto sugere, visualmente, que esta figura não seria
unidimensional, ou seja, não é uma linha dotada apenas de comprimento. Assim, a
sua dimensão estaria entre 1, dimensão da reta, e 2, dimensão do plano.
A partir da Figura 2.2, verifica-se que em cada passo, a quarta parte da curva
resultante é semelhante à curva obtida no passo anterior, existindo portanto autosemelhança nas curvas que são obtidas em nı́veis sucessivos. Portanto, no nı́vel 5 da
construção da Curva de Koch, verifica-se que a secção A0 C quatro vezes superior à
AB. Observa-se também que se houvesse uma homotetia de razão 3 na secção A0 B 0
obter-se-ia exatamente a secção AB. Este conceito de homotetia será tratado mais
claramente na seção 2.4.1. Um estudo mais detalhado de alguns fractais auto-similares
é apresentado no Apêndice A deste trabalho.
Também em 1918, Gaston Julia e Pierre Fatou, viriam a apresentar um trabalho
7
Figura 2.2: Os seis primeiros nı́veis para a construção da Curva de Koch.
8
sobre processos iterativos envolvendo números complexos, que mais tarde viriam a ser
conhecidos como “Conjuntos de Julia” [10]. Em razão da reduzida capacidade de
cálculo, então disponı́vel, não foi então possı́vel produzir qualquer imagem gráfica
destes conjuntos. No entanto, todos estes objetos matemáticos possuem algumas
caracterı́sticas comuns aos fractais. Além de conterem em si um número infinito de
elementos, apresentavam semelhança entre figuras onde uma era a ampliação da outra.
Foi, no entanto, a partir da segunda metade do século passado que questões
relativas à similitude entre uma figura e a sua ampliação começariam a aparecer
na literatura cada vez com mais freqüência. Edward Lorenz, um meteorologista
americano, dedicava-se em 1961 a previsões metereológicas. Apoiava-se em um computador, tão moderno quanto possı́vel para a época, buscando aumentar a confiabilidade
das previsões meteorológicas. Ao tentar, computacionalmente, repetir uma experiência,
se enganou nos números que deveria introduzir para o cálculo, truncando-lhe casas
decimais. Deu-se em conta que os resultados finais eram significativamente diferentes
dos que haviam sido obtidos anteriormente. Julgou, inicialmente, tratar-se de algum
outro problema que desconhecesse, mas, de fato, veio a constatar que pequenas
alterações dos valores iniciais provocavam enormes discrepâncias finais [11]. A este
fenômeno seria dado o nome de “Efeito Borboleta” pela possibilidade simbólica de o
bater de asas de uma borboleta em Pequim poder provocar um tufão em Nova Iorque.
Posteriormente, já na década de 70, James Yorke [12] veio a encontrar nos
trabalhos de Lorenz a chave para os problemas sobre os quais então se debruçava.
Cunhou assim o nome Caos. Neste surgimento de uma nova área de investigação,
contribuı́ram também os trabalhos de May e Hoppensteadt [13, 14]. Estavam assim
criadas as condições para o aparecimento de trabalhos de sı́ntese, como os de Benoit
Mandelbrot [7].
Mandelbrot, nasceu em Varsóvia em 1924 tendo se refugiado com a famı́lia em
Paris em 1936. Apesar de ter feito os seus estudos básicos de uma forma irregular,
ingressou na École Politechnique com o intuito de receber formação universitária. Em
1952 doutorou-se em Matemática pela Universidade de Paris e, em 1958, emigraria
uma vez mais. Desta feita para os Estados Unidos, iniciando uma carreira pouco
9
comum no Thomas J. Watson Research Center da IBM. Estudou a variação dos
preços de algodão, desenvolveu um trabalho relacionado com a transmissão de ruı́do
em linhas telefônicas, ensinou em Harvard e investigou a teoria dos jogos, entre outras
atividades. Em particular, Mandelbrot debruçou-se sobre um problema antigo, qual
seria o de questionar sobre qual seria efetivamente o comprimento da linha de costa
de um paı́s. Esta e outras questões estiveram na origem de toda uma teoria inovadora
que viria a culminar com a publicação do primeiro livro de Mandelbrot em 1977 [15],
que no entanto não foi bem recebido pela comunidade cientı́fica. Só em 1982, com a
publicação de “The Fractal Geometry of Nature” [7], sairia ele do anonimato.
2.2
Dimensão Fractal
Inicialmente é necessário esclarecer que ao se falar em dimensão fractal não
está se fazendo referência à dimensão euclidiana, esta usualmente denominada de
dimensão topológica. Diz-se que dois espaços topológicos têm a mesma dimensão se,
entre os pontos de um e de outro, existir uma correspondência ponto-a-ponto contı́nua
e unı́voca [16]. Pela definição de Euclides, um ponto tem dimensão 0, uma linha possui
dimensão 1, uma superfı́cie dimensão 2 e uma porção do espaço, dimensão 3. Ao se
obsevar a curva de Koch, considerando-se a geometria Euclidiana, se supõe que ela
é uma linha curva, que assim deveria ter uma dimensão topológica igual a 1. Este
tipo de curva tem um número de infinitas “dobras” que, se ampliadas, continuam
aparecendo indefinidamente e devido a este infinito detalhamento, ocupa mais espaço
que uma curva convencional (dimensão maior que 1,0), sem chegar a ocupar o espaço
do que seria uma faixa que a contém (dimensão 2). Desta maneira se faz necessário
a definição de um novo conceito de dimensão. Diferentemente do caso euclidiano,
o termo dimensão fractal, refere-se à dimensão espacial sendo portanto associado às
irregularidades de um dado objeto no espaço que o contém.
O termo fractal aplica-se, em geral, a construções muito diversas, tanto nas
formas ditas como abstratas, como nas formas inerentes à natureza e que são objeto de
estudo da Fı́sica. Em razão disto, existe a necessidade de se estabelecer uma definição
10
do que é um fractal e que atenda à todas as espécies de construção. Para Mandelbrot
[7], “Um dado conjunto A constitui um fractal se, em A, D > Dt , sendo D a dimensão
Fractal e Dt a dimensão topológica do conjunto A”. Kenneth Falconer [16], propõe
uma definição menos rigorosa em termos das caracterı́sticas dos fractais. Uma destas
caracterı́sticas seria a complexidade infinita, ou seja, o fato que sucessivas ampliações
de um fractal levam, indefinidamente, a mais detalhes. Para figuras geométricas
convencionais, como por exemplo a circunferência, isto não ocorre. Uma estrutura
fractal não pode ser descrita de maneira simples, por exemplo, uma vez que o fractal é
construı́do por processos iterativos. Além disso, os fractais possuem auto-similaridade,
o que consiste em se poder obter réplicas de menores porções de uma figura através
de sua ampliação.
Para definir, de maneira generalizada, a dimensão fractal de um conjunto,
considera-se o espaço métrico Rn . Um espaço métrico é um par (M, d), onde M é um
conjunto não vazio e d(métrica) uma função d : M X M → R em que, para todos os
pontos x e y de M são satisfeitas as relações:
(i) d(x, y) ≥ 0, (positiva)
(ii) d(x, y) = 0 ⇔ x = y, (não degenerada)
(iii) d(x, y) = d(y, x), (simétrica)
(iv) d(x, z) ≤ d(x, y) + d(y, z), (desigualde triangular)
Na reta definimos um intervalo como um segmento de reta, enquanto que no
Rn o intervalo é definido como uma esfera de raio ε, centrada em x. Esta esfera será
representada por B(x, ε), sendo então uma bola aberta definida como : B(x, ε) = {y ∈
Rn /d(x, y) < ε}, onde d(x, y) é a métrica do espaço Rn . Já uma bola fechada é
definida como: B(x, ε) = {y ∈ Rn /d(x, y) ≤ ε}. Representa-se o menor número de
bolas fechadas de raio ε necessário para ter uma cobertura A, com A ⊂ Rn e ε > 0,
por N (A, ε).
Um conjunto A, terá dimensão D se [4]:
11
N (A, ε) ∼
= Cε−D
(2.5)
sendo C uma constante positiva.
Daı́ segue-se que:
µ ¶D
N (A, ε) ∼ 1
=
C
ε
(2.6)
Calculando-se o logaritmo neperiano dos dois membros da relação 2.6, tem-se:
µ ¶D
N (A, ε) ∼
1
ln(
) = ln
C
ε
(2.7)
E, portanto:
ln(N (A, ε) − ln(C)
D∼
.
=
ln( 1ε )
Como
ln(C)
ln( 1ε )
(2.8)
tende a 0 à medida em que ε → 0, segue-se a seguinte definição:
Definição 1 Seja dado um conjunto compacto [17] A ⊂ Rn . Para cada ε > 0 ,
N (A, ε) denota o menor número de bolas fechadas de raio ε necessário para cobrir A.
Então, a dimensão fractal deste conjunto é definida por:
"
ln(N (A, ε))
D = lim
ε→0
ln( 1ε )
#
(2.9)
Caso a variável contı́nua ε seja substituı́da por uma variável discreta, tem-se o
seguinte teorema:
Teorema 1 Seja A ⊂ Rn . Dado εn = Crn em R, onde 0 < r < 1, C > 0 e o inteiro
n={1,2,3,...}, então A, tem dimensão fractal será dada por:
"
ln(N (A, εn ))
D = n→∞
lim
ln( ε1n )
#
(2.10)
Prova 1 Sejam r, C reais e a seqüência E = {εn : n = 1, 2, 3, ...}, conforme descritos
no teorema. Definindo f (ε) = M ax{εn ∈ E; εn ≤ ε}. Assumindo que ε ≤ r, então:
12
f (ε) ≤ ε ≤
f (ε)
r
(2.11)
Determinando-se os inversos de cada termo da desigualdade e aplicando-se na
seqüência o logaritmo neperiano, ter-se-á:
1
1
r
≥ ≥
f (ε)
ε
f (ε)
(2.12)
e
Ã
1
ln
f (ε)
!
µ ¶
r
1
≥ ln(
)
ε
f (ε)
≥ ln
(2.13)
o que levará à seguinte desigualdade:
ln[N (A, f (ε))] ≥ ln[N (A, ε)] ≥ ln[N (A,
f (ε)
)]
r
(2.14)
Como lnx é uma função crescente positiva para x ≥ 1, e uma vez que:
1
1
)
ln( f (ε)
≤
1
ln
³ ´ ≤
1
ε
1
(2.15)
r
ln( f (ε)
)
segue-se que:
ln[N (A, f (ε)
)]
ln[N (A, f (ε))]
ln[N (A, ε)]
r
³ ´
³
´
≥
≥
r
1
1
ln( f (ε) )
ln ε
ln f (ε)
(2.16)
Calculando o limite do lado esquerdo da desigualdade 2.16 com N (A, ε) → ∞
e ε → 0, tem-se:
lim {
n→∞
ln[N (A, f (ε))]
ln[N (A, εn )]
} = n→∞
lim {
}
r
ln( f (ε) )
ln( εrn )
(2.17)
o que levará a:
lim {
n→∞
ln[N (A, εn )]
ln(r) + ln
³
1
εn
´ } = lim {
n→∞
pois ln(r) é constante.
13
ln[N (A, εn )]
ln
³
1
εn
´
}
(2.18)
Calculando agora, o limite do lado direito da desigualdade 2.16 tem-se:
ln[N (A, f (ε)
)]
ln[N (A, εn−1 )]
r
lim
{
} = limn→∞ {
}
1
n→∞
ln( f (ε) )
ln( ε1n )
(2.19)
e portanto:
lim {
ln[N (A, εn−1 )]
n→∞
ln( 1r )
+ ln
³
1
εn−1
´ } = lim {
n→∞
ln[N (A, εn−1 )]
ln
³
1
´
} = lim {
n→∞
εn−1
ln[N (A, εn )]
ln
³
1
εn
´
}
(2.20)
Como ε → 0 por ambos os lados da desigualdade, aproximando-se do mesmo
valor, recorrendo-se ao teorema do confronto [18] do cálculo, o limite do termo do
meio existe e é igual ao limite dos extremos, logo:
"
ln(N (A, εn ))
D = lim
n→∞
ln( ε1n )
2.3
#
(2.21)
Fractais Auto-Similares e Fractais Auto-Afins.
As estruturas com geometria fractal apresentam uma auto-similaridade que
pode ser exata ou estatı́stica, sendo então denominadas de fractais determinı́sticos ou
randômicos, respectivamente. Por auto-similaridade exata entende-se a invariância
da estrutura após uma transformação isotrópica. Se tomarmos um objeto S, formado
por um conjunto de pontos R = {x1 , x2 , x3 , ...}, a tranformação similar com um
fator de escala b, muda as coordenadas dos pontos para bR = {bx1 , bx2 , bx3 , ...}.
Logo, o conjunto S, formado pelos pontos de coordenadas R, é auto-similar, se este
resulta invariante após a referida transformação. Para exemplificar, considere-se um
fractal denominado Triângulo de Sierpinski, como mostrado na Figura 2.3. Este
fractal é aproximadamente 40 anos “mais jovem” que o conjunto de Cantor e foi
apresentado pelo matemático polonês Waclaw Sierpinski [19]. Sua construção básica
começa com um triângulo equilátero, totalmente preenchido. Inicialmente toma-se os
pontos médios dos três lados que, juntos com os vértices do triângulo original, formam
quatro triângulos congruentes. Em seguida, retira-se o triângulo central, concluindo14
se assim o processo básico de construção. Tem-se então três triângulos congruentes,
cujos lados medem metade do lado do triângulo original. Repete-se, com cada um
destes três triângulos, o processo descrito, tantas vezes quanto se desejar. Desta
maneira começa-se com um triângulo e geram-se sequencialmente 3, 9, 27, 81, 243,
... triângulos, segundo uma mesma lei de formação. Sua geometria triangular é
constituı́da por triângulos menores que são cópias perfeitas da forma completa da
figura. Assim, ao ampliar-se (“zoom”) uma parte qualquer, se terá algo idêntico à
figura como um todo. Neste caso, esta é uma figura fractal auto-similar exata, ou
simplesmente denominada auto-similar.
Figura 2.3: O triângulo ADE, com todo seu conteúdo, é uma redução exata do
triângulo ABC. O mesmo se pode dizer com relação aos triângulos CDF e de BEF.
Existem contudo, fractais que são igualmente formados por mini-cópias, mas
estas são anisotrópicas, ou seja, não mantêm fixas as proporções originais. Ao se passar
de uma escala para outra maior, observa-se que o tamanho destas cópias não aumenta
uniformemente em todas as direções espaciais. Neste caso, os fractais são chamados
de auto-afins. Na Figura 2.4, tem-se um fractal auto-afim gerado a partir da diagonal
15
de um quadrado que é dividido em quatro partes iguais horizontalmente e reproduzida
de acordo com a Figura 2.4. No próximo nı́vel, repete-se o mesmo procedimento com
os quatro segmentos dispostos e assim sucessivamente. Por outro lado, diz-se que
esta auto-afinidade é definida por suas caracterı́sticas estatı́sticas, ou seja, a fronteira
da figura mantém uma correlação estatı́stica quando observada em diferentes escalas.
Portanto, para se definir auto-afinidade, utiliza-se uma transformação anisotrópica,
ou seja, diferentes fatores de escala são associados a diferentes direções espaciais
para obtenção de uma configuração estatisticamente semelhante. Sendo assim, bR =
{b1 x1 , b2 x2 , b3 x3 , ...}.
Figura 2.4: Exemplo de um fractal auto-afim. Neste caso o fator de escala para as
abscissas é 4 e para as ordenadas é 2.
16
2.4
Métodos para Determinação da Dimensão Fractal
2.4.1
Homotetia e Dimensão Fractal
Intuitivamente, a definição de dimensão fractal se torna simples adotando-se
como ponto de partida o conceito de homotetia [20]. Homotetia significa ampliação
ou redução de qualquer ente geométrico, podendo ser figuras planas, como triângulos,
quadriláteros, cı́rculos, etc, ou espaciais como cubos, pirâmides, esferas, ou figuras
com geometria irregular.
Para definir-se homotetia necessita-se da escolha de um centro de homotetia,
O, e da razão de homotetia ,r, sendo que esta pode ser qualquer número real. Um
novo ponto, X 0 , em relação a O, é gerado através do ponto X inicial, por homotetia,
através da equação:
−−→0
−−→
OX = r.OX
(2.22)
Esta definição pode ser visualizada pelo exame da Figura 2.5. Este tratamento
vetorial na definição de homotetia serve para garantir que dois pontos homotéticos
entre si estão alinhados com o centro de homotetia. Além disso, se r > 0, X e X 0
estão em uma mesma semi-reta onde O é um extremo. Se r < 0, então X e X 0 estão
em semi-retas distintas, sendo que O é extremo das duas.
Tomando-se, por exemplo, um segmento de reta, podemos dividı́-lo em p partes
iguais, semelhantes ao segmento original, porém reduzidas em uma certa razão r.
Obviamente o número n de segmentos obtidos relaciona-se com a razão de semelhança,
r, de modo que n = 1r . Para um quadrado, fazendo-se uma divisão de cada um dos
seus lados em p partes iguais, obtem-se p2 quadrados semelhantes ao original, de modo
que p2 = n e, portanto, n =
1
.
r2
Finalmente, para um cubo, dividindo-se cada uma
de suas arestas em p partes iguais, a razão de semelhança entre os lados obtidos e
o original é de
1
p
e então n =
1
.
r3
Assim, considerando-se que a dimensão espacial
das figuras tradicionais é igual à sua dimensão topológica, um segmento de reta terá
17
Figura 2.5: Figura representativa da operação de homotetia. O ponto C é gerado
através de OC’=r.OC, com r < 1,0. Tem-se A’C’//AC e B’C’//BC. Então os
triângulos, 4(A0 B 0 C 0 ) e 4(ABC), são semelhantes.
dimensão 1, um quadrado terá dimensão 2 e um cubo dimensão 3. Deste modo, é
razoável afirmar que:
n=
1
rD
(2.23)
onde D é a dimensão espacial, r é razão de semelhança e n o número de réplicas da
figura.
Para os fractais auto-similares, ou seja, aqueles cuja forma geométrica independe do comprimento de escala, basta analisar o que ocorre de uma etapa para outra da
construção, uma vez que isto reflete o que acontece em todo conjunto. Calculando-se
o logaritmo neperiano em ambos os membros da equação 2.23, a dimensão tem a
forma:
D=−
ln n
ln r
(2.24)
onde n representa o número de partes restantes após uma etapa, ou geração da
construção do fractal, sendo r a razão de semelhança destas partes com a configuração
original.
Tomando-se, como exemplo, o Conjunto de Cantor, para cada geração tem-se
dois segmentos, que serão posteriormente seccionados em três partes, logo n = 2. Já a
18
razão de semelhança desses segmentos com a figura original é 13 . Portanto, a dimensão
Fractal do conjunto de Cantor é dada por:
D=−
2.4.2
ln 2
ln 2
≈ 0, 63
1 =
ln 3
ln( 3 )
(2.25)
Método de Contagem de Caixas.
Caso o fractal não apresente auto-similaridade exata, i.e., tem-se um fractal
auto-afim, o cálculo da correspondente dimensão fractal torna-se mais difı́cil, ou
muitas vezes impossı́vel, por uma metodologia baseada na homotetia. Para uma
situação como esta, um método que antes era freqüentemente utilizado é o de Contagem de Caixas [21]. Atualmente, sabe-se que este método tem boa utilidade para
casos em que o fractal seja auto-similar [22]. Este método tem sido utilizado para
o cálculo da dimensão fractal de nuvens, de formações em áreas fotografadas por
satélites, entre outras [23]. O método consiste em “cobrir” a figura, da qual se deseja
determinar a dimensão fractal, com uma malha de quadrı́culas de lado r, conforme
mostrado na Figura 2.6. Sendo n o número mı́nimo de quadrı́culas, de lado r, que
contêm ao menos 1 ponto da figura, e δ o lado da moldura que escolhe-se para nela
inserir a figura, então:
Ã !D
n=
δ
r
(2.26)
Como exemplo da aplicação do método, considere-se que a figura da qual
se quer determinar a dimensão seja um quadrado.
Dividindo-se o mesmo em 4
quadrı́culas iguais, cada uma delas terá a metade do lado que tinha o quadrado
anterior. Logo p = 2 e
δ
r
= 2. Então, verifica-se claramente que D = 2 o que condiz
com a dimensão de um quadrado cheio. Da mesma forma, dividindo-se o mesmo em
16 quadrı́culas iguais, cada uma delas terá um quarto do lado que tinha o quadrado
anterior. Logo p = 4 e
δ
r
= 4. Do mesmo modo, verifica-se que D = 2.
Para que a determinação da Dimensão Fractal seja mais precisa, é necessário
que a malha seja tal que o lado r das quadrı́culas seja tão pequeno quanto se queira,
19
de modo que, no limite para r tendendo a 0, tem-se:
D = lim
r→0
ln(n)
ln( rδ )
(2.27)
É possı́vel implementar o método de contagem de caixas computacionalmente,
através de algoritmos simples em diversas linguagens de programação. Primeiramente
utiliza-se uma grade de quadrı́culas de lado r1 . Conta-se o número de quadrı́culas que
contenham pelo menos um ponto da figura. O processo é repetido para uma malha
de quadrı́culas de lado r2 , sendo r2 < r1 , contando-se o número n2 de quadrı́culas que
contenham ao menos um pixel da figura.
Utilizando-se a relação 2.26 para para n1 quadrı́culas de lado r1 , então:
Ã
n1 =
δ
r1
!D
(2.28)
Para n2 quadrı́culas de lado r2 , tem-se analogamente:
Ã
n2 =
δ
r2
!D
(2.29)
Dividindo-se membro a membro as relações 2.28 e 2.29, e isolando-se o valor
da dimensão fractal, ter-se-á:
D=
ln n1 − ln n2
ln r2 − ln r1
(2.30)
Assim, torna-se possı́vel calcular a dimensão fractal, pois os valores de n1 , n2, r1
e r2 são conhecidos. O método pode se tornar mais preciso construindo-se m grades
diferentes, i.e., com quadrı́culas de tamanhos diferentes que cobrem a figura. Cada
quadrı́cula possui lado ri , sendo ni o número de quadrı́culas da i-ésima grade que
contêm pelo menos um ponto da figura (i = 1, 2, 3, ..., m).
Para cada uma dessas grades, conta-se quantas quadrı́culas contêm ao menos
um ponto do objeto e, a partir disso, constróe-se um gráfico do número de quadrı́culas
contendo partes do objeto em função do tamanho das quadrı́culas. O valor absoluto
20
Figura 2.6: A imagem de uma galáxia sendo utilizada para o procedimento de
contagem de caixas. As imagens mostram o processo sucessivo de contagem
para caixas com diferentes tamanhos 75x75, 40x40, 22x22, 12x12 e 3x3
pixels. As caixas brancas contêm o contorno e as pretas, não interceptam o
mesmo.(http://www.ees.nmt.edu/ davew/P362/boxcnt.htm).
Figura 2.7: A inclinação da reta ajustada determina a dimensão fractal da imagem
mostrada na Figura 2.6, neste caso 1, 54 ± 0, 04.
21
da inclinação da reta assim definida em escala logarı́tmica dá uma estimativa da
dimensão fractal do objeto investigado. A reta a + b log d representada na Figura 2.7
para o caso mostrado na Figura 2.6, que melhor se ajusta a esse conjunto de pontos,
a partir do método dos mı́nimos quadrados, tem por parâmetros:
µm
P
a=
i=1
¶ µm
P
log ni .
¶
2
(log ri )
i=1
m.
µm
P
−
µm
P
i=1
log ri
¶
(log ri
i=1
)2
−
µm
P
¶µ m
P
¶
(log ri )(log ni )
i=1
(2.31)
¶2
log ri
i=1
e
m.
b=
µm
P
¶
(log ri )(log ni ) −
i=1
m.
2.4.3
µm
P
¶
(log ri
i=1
)2
µm
P
i=1
−
log ri
µm
P
i=1
¶µ m
P
i=1
¶2
¶
log ni
(2.32)
log ri
Método RMS
Este método é particularmente simples [24], pelo menos conceitualmente,
permitindo identificar se perfis/superfı́cies são fractais auto-afins. Tome-se como
exemplo o caso de um perfil. Calcula-se o desvio padrão dos valores das alturas,
em relação à altura média, dentro de uma janela quadrada de largura w, como
mostrado na Figura 2.8.
A precisão do método é melhorada tomando-se várias
janelas. Todas de mesma largura, ou de mesma área para uma superfı́cie, dispostas
em diferentes posições do perfil. A maneira de se distribuir estas janelas pode ser
sistemática ou randômica. No primeiro caso, as sucessivas janelas são construidas com
posicionamento da origem seqüencialmente bem definido e, no segundo, esta origem é
escolhida aleatoriamente. Calcula-se o desvio padrão correspondente a cada janela e,
a seguir, o valor médio dos desvio padrão para as múltiplas janelas. O procedimento
é então repetido para conjuntos de janelas com tamanhos diferentes e o desvio padrão
médio é verificado em função do tamanho da janela. O RMS (Root Mean Square) ou
rugosidade é freqüentemente utilizado para caracterizar morfologicamente superfı́cies,
sendo definido como:
22
v
nw u
1 X
u 1 X
t
(zj − < z >)2
RM S(w) =
nw i=1 mi j∈w
(2.33)
Na equação acima, nw é o número total de janelas de comprimento w, mi é
o número de pontos dentro da i-ésima janela, zj corresponde à altura de cada ponto
de uma mesma janela e < z > é a altura média deste conjunto de dados contidos na
i-ésima janela, cujo comprimento é w.
Figura 2.8: Procedimento para o cálculo do RMS com janelas de lado w, sobre perfis
topográficos. Sucessivas janelas de um mesmo tamanho são construı́das, podendo a
sua distribuição ao longo do perfil ser sistemática ou aleatória. (J.G.V. Miranda,
Análisis Fractal Del Microrrelieve del Suelo. La Coruña: Espanha, 2000. p. 113 Tese
de Doutorado)
No caso de duas dimensões teremos que a janela que abrange parte de um
perfil é substituı́da por um quadrado de lado w. No caso de perfis ou superfı́cies com
geometria fractal, o RMS mudará dependendo do comprimento de escala w utilizado
no cálculo. Isto se deve ao fato do perfil/superfı́cie ser invariante após uma mudança
de escala. Para geometrias auto-afins [4], a rugosidade varia com o comprimento de
escala utilizado de acordo com a relação:
RM S(w) ≈ wα
23
(2.34)
Na expressão acima, calculando-se o logaritmo de ambos os termos obtém-se
a linearização em que o valor do expoente de rugosidade, α, é expresso em função
da inclinação da reta. Esta reta é ajustada entre duas escalas, uma correspondente à
escala inicial e a outra à escala final.
É importante discutir também, como se calcular a dimensão fractal a partir
do cálculo do expoente de rugosidade já determinado pelo método RMS. Para a
determinação da dimensão fractal de uma superfı́cie auto-afim, parte-se da definição
do que se entende com uma superfı́cie assim caracterizada. Seja uma variável z,
dependente de x e y. Considere-se um fator de escala, bα , aplicado a z, distinto do
fator aplicado a x e y, nestes casos o fator b. Esta relação pode ser escrita de acordo
com:
z(bx, by) ∼ bα z(x, y)
(2.35)
A relação acima mostra o fato de que uma superfı́cie auto-afim pode ser
escalada com diferentes fatores nas direções x, y e z. Isto quer dizer que, no caso, se
reescala as direções x e y por um fator multiplicativo b, (x → bx e y → by), ao tempo
em que a direção z foi reescalada por um fator multiplicativo bα (z → bα z). Para o
caso particular em que α = 1 a transformação seria auto-similar.
Uma importante conseqüência da equação 2.35 é o comportamento de lei de
potência entre a diferença dos valores das alturas z(x, y), nos pontos (x1 , y1 ) e (x2 , y1 ),
e a distância entre estes últimos. No caso de duas dimensões, esta relação é tal que,
definindo-se δ ≡ |z(x1 , y1 ) − z(x2 , y1 )| e l ≡ |(x1 − x2 )|, tem-se [4]:
δ ∼ lα
(2.36)
Consideremos, coerentemente com a equação anterior, um perfil, para um valor
constante de y (y1 = constante), de uma superfı́cie auto-afim, definida de modo que
os valores de x e de y pertençam ao intervalo [0,1]. Primeiramente, ao dividir-se
o domı́nio da direção x em N segmentos iguais, o comprimento de cada segmento
unitário é dado por l =
1
.
N
No intervalo horizontal de tamanho l a altura muda
24
de acordo com a equação 2.36, de modo que se requerem
δ
l
∼ lα−1 caixas de lado l
para cobrir o perfil auto-afim. Assim, para N segmentos, o número total de caixas
requeridas para recobrir o perfil, N (l), será dada por:
N (l) ∼ lα−2
(2.37)
Como para fractais em geral tem-se que N (l) ∼ l−Df , então a dimensão fractal
para um perfil auto-afim será:
Df = 2 − α
(2.38)
Para uma superfı́cie auto-afim:
Df = 3 − α
2.4.4
(2.39)
Método Semivariograma
O semivariograma é uma função que relaciona a semivariância com o valor
dos incrementos de distância. Serve para descrever a variabilidade espacial de um
fenômeno. Existem uma série de modelos teóricos que podem descrever uma função
monotonicamente crescente assim como o semivariograma, que pretende sintetizar
toda a informação que se conhece a respeito da correlação espacial da variável de
interesse. Este método é baseado no cálculo da dependência espacial mediante estimativas da semivariância nas diferentes escalas de medida. A semivariânçia em função da
distância proporciona um semivariograma, termo que também se utiliza para designar
tal método. Por definição a equação do semivariograma para uma função contı́nua
Z(x) é expressa por:
1
γ ∗ (x) = h[Z(x) − Z(x + h)]2 i
2
(2.40)
Na equação 2.40, h é a distância a partir da origem ou escala e o sı́mbolo h i
significa uma integral sobre todo o domı́nio da função Z. Para o caso deste trabalho
25
a distribuição de pontos para a superfı́cie considerada é discreta e assim, definindose {X}, {Y } e {Z} os conjuntos que contém as coordenadas Xi e Yi que permitem
localizar cada ponto no espaço bi-dimensional e Zi a altura do mesmo, pode-se definir
uma função discreta Z que relacione os elementos dos conjuntos {X} e {Y } com os do
conjunto {Z}.Denotando tal função por Z(xi , yi ), a semivariância pode ser definida
por:
X
1 n(h)
γ(h) =
[Z(xi , yi ) − Zh ]2
2n(h) i=1
(2.41)
Na equação 2.41, h é a escala ou distância média, calculada no intervalo [h−∆h,
h + ∆h], sendo 2∆h a distância de separação entre as escalas no semivariograma; o
termo [Z(xi , yi ) − Zh ] representa a diferença de altura entre pontos separados pela
distância média h e n(h) é o número de pontos considerados em cada intervalo.
Para mostrar a dependência da semivariância com o expoente de rugosidade,
analisemos o caso unidimensional de uma partı́cula com movimento browniano que se
move ao longo de uma linha dando saltos de comprimento +ε e −ε sendo cada salto
realizado num intervalo de tempo denotado por τ . Definamos o intervalo de tempo τ
fixo e o comprimento do passo ε dado por uma distribuição normal:
³
1
e
P (ε, τ ) = √
4πF τ
2
´
ε
− 4F
τ
(2.42)
onde o parâmetro F é o coeficiente de difusão.
Os comprimentos ε são definidos de forma que uma seqüência de passos {ε}
constituem um conjunto de variáveis independentes, aleatórias e Gaussianas. Logo,
a probabilidade de encontrar um comprimento de passo ε em um intervalo entre ε e
ε + dε é dada por P (ε, τ )dε e a variância do processo está definida por:
+∞
Z
2
ε2 P (ε, τ )dε = 2F τ
hε i =
(2.43)
−∞
Naturalmente, ao se somar uma seqüência de passos se obtém a distância da
partı́cula à origem de modo que:
26
X(t = nτ ) =
n
X
εi
(2.44)
i=1
Considerando-se dois incrementos ε0 e ε00 estatisticamente independentes a
distribuição de probabilidades total será o produto das densidades de probabilidade e
portanto:
P (ε0 , ε00 , τ ) = P (ε0 , τ )P (ε00 , τ )
(2.45)
Integrando-se sobre todas as seqüências de combinações possı́veis de ε0 e ε00
obtém-se a densidade de probabilidade para o incremento ε = ε0 + ε00 de acordo com
a expressão:
³
+∞
Z
1
P (ε − ε , τ )P (ε , τ )dε = √
e
4πF 2τ
0
P (ε, 2τ ) =
−∞
0
0
2
− 4Fε 2τ
´
(2.46)
Observa-se portanto que a distribuição 2.46 possui média nula e variância
hε2 i = 4F τ, ou seja, incrementada por um fator multiplicativo de 2 em relação à
expressão 2.43. Tal situação pode ser generalizada para intervalos de tempo bτ de
modo que:
³
1
P (ε, bτ ) = √
e
4πF bτ
2
− 4Fε bτ
´
(2.47)
A partir da expressão 2.47, naturalmente hε2 i = 2F bτ de modo que convenientemente h[X(t+bτ )−X(t)]2 i = 2F bτ . Isto significa que a variância depende da escala
de observação. Denotando-se h = bτ , como o fator de escala utilizado, evidencia-se
uma dependência da função semivariância para o movimento browniano com a escala
utilizada de modo que:
γ(h) ∝ h
(2.48)
Para um movimento browniano fracionário generalizado a relação h[X(t+bτ )−
X(t)]2 i = 2F bτ resulta em:
27
h[Xα (t + bτ ) − Xα (t)]2 i = 2F h2α
(2.49)
onde α é um parâmetro que pode variar entre 0 < α < 1. Note que para um
movimento browniano simples com os incrementos estatisticamente independentes,
α = 0, 5 por 2.48, de modo que a variável posição pode ser denotada por X0,5 (t) =
X(t).
2.4.5
Análise de Fourier
A análise de Fourier ou análise espectral de uma série temporal, consiste em
expressar uma função f (t), considerando os fenômenos investigados, de preferência
periódicos, como uma somatória de funções trigonométricas de senos e cossenos. As
séries em princı́pio são infinitas, mas na prática, quando se trabalha com espaços
discretos, a freqüência mais elevada com importância, corresponde ao espaçamento
dos “pixels” (amostras unitárias que compõem uma imagem - picture elements) na
imagem ou ao longo da elevação de um perfil. Do mesmo modo, a menor freqüência
de interesse é ajustada pela largura da imagem ou o comprimento de um perfil. Para
uma função do tempo a equação de Fourier se torna:
f (t) =
kX
max
ak .sen(2πkt − ϑk )
(2.50)
kmin
Substituindo t por x, a expressão 2.50 pode ser utilizada para perfis com
elevações. Os coeficientes ak dão a magnitude da k − ésima freqüência, e os valores ϑ
informam as fases de cada termo da série. Uma caracterı́stica muito útil do método de
Fourier é que os coeficientes para cada termo não mudam quando outros termos são
adicionados às séries. Porém, como uma seqüência finita de freqüências são usadas, o
resultado da soma é uma “melhor aproximação”. Comumente, o comportamento do
quadrado da amplitude como função da freqüência resulta no espectro de potências
do perfil original.
Como um perfil fractal, por definição, inclui informações de todas as freqüências,
um pensamento preliminar concluiria que este método seria difı́cil para se aplicar
28
a fractais. Contudo o comportamento verificado é surpreendentemente simples. O
espectro de amplitude mostra uma variação linear entre o logaritmo do quadrado da
amplitude e o logaritmo da freqüência. A inclinação desta reta está relacionada com a
dimensão fractal do perfil, de modo que esta técnica se torna aplicável à dados autoafins ou auto-similares. Sendo β a inclinação, no ajuste linear, a dimensão fractal
para perfis será dada por [25]:
Df =
(4 + β)
2
(2.51)
A análise de Fourier também pode ser executada para superfı́cies. Neste caso, a
transformação consiste dos dados de amplitude e fase, convertidos numa disposição bidimensional. Quando a escala da imagem de uma superfı́cie fractal é a transformada
de Fourier, o resultado esperado é que o logaritmo do quadrado da amplitude cairá
linearmente com o logaritmo da freqüência [25]. Para medir finalmente a dimensão
fractal para uma superfı́cie fractal, tem-se:
Df =
2.5
(6 + β)
2
(2.52)
Simulação Computacional de Superfı́cies Fractais
Nesta seção discutiremos duas metodologias distintas para a geração, mediante simulação computacional, de superfı́cies com geometria fractal. Superfı́cies assim geradas
serão utilizadas no desenvolvimento deste trabalho, funcionando como fronteira de
uma região onde se deseja determinar o comportamento do campo e do potencial
elétrico.
2.5.1
Técnica FBM
Um modelo matemático para a geração de um objeto auto-afim (também
chamado de estatisticamente auto-similar), é o que resulta da técnica denominada de
29
FBM (Fractional Brownian Motion), técnica esta que apresenta propriedades tı́picas
daquelas caracterı́sticas do Movimento Browniano [26, 24]. Em uma dimensão, o FBM
é um processo randômico Z(t), capaz de gerar uma variável Z como função aleatória
de um parâmetro t. A variável Z apresenta incrementos Gaussianos, Z(t2 ) − Z(t1 ),
onde a variância destes é proporcional a |t2 − t1 |2α com 0 < α < 1, ou seja, os
incrementos de Z são auto-afins, conforme sub-seção 2.4.4. Isto significa que tomandose t0 = 0 e Z(t0 ) = 0, as duas funções randômicas Z(rt) e rα Z(t), são estatisticamente
indistinguı́veis, com r sendo um fator de escala. Para um dado número Z0 , tem-se que
os pontos t que satisfazem Z(t) = Z0 , constituirão pontos pertencentes a um fractal
auto-afim e com dimensão fractal Df = 1 − α (ver Figura 2.9).
Figura 2.9: Esquema mostrando a interseção da reta Z0 com os pontos Z(t) pertencentes a um
movimento browniano fracionário. Este conjunto de pontos de interseção constitui o que se chama
de Seção de Poincaré, com dimensão fractal Df = 1 − α.
A generalização da técnica FBM para dimensões superiores é um processo
multidimensional Z(t1 , t2 , ..., tn ), onde Z, tem incrementos estacionários e anisotrópicos, isto é, todos os pontos (t1 , t2 , ..., tn ) e todas as direções são estatisticamente
equivalentes.
Alguns algoritmos tem sido desenvolvidos para gerar computacionalmente fractais randômicos e compará-los com estruturas naturais que apresentam geometria
fractal. Estes métodos têm tido a sua qualidade julgada mediante a comparação do
que é gerado computacionalmente vis-a-vis o que se tem na natureza, bem como com
respeito ao seu custo computacional. Muito recentemente, com a disponibilidade de
30
computadores com maior rapidez de processamento e de armazenamento, bem como
de programas gráficos, o controle de propriedades locais e gerais destes tais fractais
computacionalmente gerados se tornou possı́vel. Tem-se como exemplo o controle local
da dimensão fractal que é utilizado para modelar os vales cercados por montanhas
ásperas em uma cena de paisagem.
Um dos primeiros e mais conhecidos métodos de geração de fractais FBM é o
“Midpoint Displacement Method” [27]. Nesta metodologia, assume-se que os valores
Z(0) = 0 e Z(1) sejam dados, com Z(1) podendo ser obtido num processo randômico
cuja variância é σ 2 . O intervalo [0, 1] é particionado em dois sub-intervalos [0, 12 ], [ 12 , 1]
e Z( 21 ) é definido como a média entre Z(0) e Z(1), adicionado de um deslocamento
D1 , isto é:
1
1
Z( ) = (Z(0) + Z(1)) + D1
2
2
(2.53)
O deslocamento D1 é computado como uma variável randômica Gaussiana, com
variância ∆21 que é proporcional a
σ2
.
22α
Esta propriedade, passı́vel de demonstração,
contém o expoente de rugosidade, α, o que permite gerar uma superfı́cie (que em uma
dimensão corresponde a um perfil), com o expoente de rugosidade que se deseja. O
processo é repetido para cada um dos dois intervalos restantes, de tal forma que:
1
1
1
Z( ) = (Z(0) + Z( )) + D2,1
4
2
2
(2.54)
1
1
3
Z( ) = (Z( ) + Z(1)) + D2,2
4
2
2
(2.55)
Os deslocamentos D2,1 e D2,2 acima, são Gaussianos com variância ∆22 proporcional a
σ2
,
(22 )2α
sendo que os dois valores de D2 nas equações acima, podem ser
diferentes. O processo iterativo, para cada novo intervalo, continua, com deslocamentos Dn,i (i = 1,...,2n) e variâncias ∆2n proporcionais a
σ2
(2n )2α
para o n-ésimo estágio.
Neste trabalho, tratou-se o caso em que Z(t1 , t2 ) é considerado como uma altura definida sobre o ponto (t1 , t2 ) de um plano. O resultado é uma superfı́cie fractal
auto-afim, com dimensão fractal Df = 3 − α. Vale observar que no processo iterativo
31
de geração o fator α é mantido constante podendo ser previamente definido.
Os conjuntos de pontos {(t1 , t2 ) satisfazendo Z(t1 , t2 ) = Z0 } são coleções de
curvas estatisticamente auto-semelhantes, com dimensões fractais Df = 2 − α (ver
Figura 2.10).
Figura 2.10: Figura mostrando os pontos (t1 , t2 ) que satisfazem Z(t1 , t2 ) = Z0 . Este conjunto de
pontos é chamado de Ilhas de Korcak e representa a mesma idéia da seção de Poincaré, extendida
contudo para o caso de superfı́cies. A dimensão fractal deste conjunto de pontos de interseção possui
dimensão fractal Df = 2 − α.
Um exemplo simples de um método que resulta em uma superfı́cie cujas alturas
são geradas randomicamente, Z(t1 , t2 ) é o que corresponde a ter-se, inicialmente, um
triângulo equilátero. Este é subdividido em quatro pequenos triângulos, também
equiláteros (Figura 2.11).
Figura 2.11:
Geração de um fractal usando o método “midpoint displacement”.
As
alturas iniciais xA , xB e xC dos vértices do triângulo original e os deslocamentos dos vértices
dos triângulos descendentes, yA , yB e yC , são mostrados para os dois primeiros nı́veis.
(http://citeseer.ist.psu.edu/470996.html)
Os vértices recém criados, associados às metades dos lados do triângulo original,
são deslocados verticalmente por valores randômicos, como no caso unidimensional.
32
Um processo similar é repetido para cada um dos triângulos menores e assim para
cada um dos seus “descendentes”, até que a interação que se deseja seja atingida. O
número de pontos para os quais é definida uma altura depende do número de iterações
que se escolheu.
2.5.2
Deposição Balı́stica
A morfologia e a evolução de uma superfı́cie, no crescimento de filmes finos, tem
sido objeto de interesse para pesquisas de cunho teórico e experimental. Um modelo
particularmente simples é o de deposição balı́stica, que foi originalmente formulado
como um modelo de sedimentação e tem sido extensivamente estudado como um
modelo de crescimento de filmes finos a baixas temperaturas [28]. O algoritmo de
crescimento para a superfı́cie balı́stica, considera, no processo de deposição, correlações
laterais e perpendiculares, de tal forma que uma dada partı́cula, depositada sobre um
substrato inicialmente descorrelacionado, atinge uma altura que poderá ser igual ou
maior à da correspondente imediatamente vizinha.
Figura 2.12: Esquema representando a formação de um perfil balı́stico. Neste modelo as partı́culas
fixam-se onde atingem a superfı́cie por contato lateral ou perpendicular, sem qualquer difusão. Na
Figura, h e i, representam respectivamente os valores de altura de cada sı́tio e do comprimento
correspondente, de modo que o tamanho total do sistema se dá para i = L.
No instante t = 0, a interface é plana, ou seja, a altura de cada sı́tio é
considerada nula. Matematicamente, h(i, t = 0) = 0 para i = 1, 2, ..., L. Em algum
33
instante t, por escolha randômica de um ponto i da rede, faz-se crescer h(i, t) para
h(i, t + 1). Isto é feito de tal modo que h(i, t + 1) = max[h(i − 1, t), h(i, t) + 1, h(i +
1, t)]. Para um sistema como este, uma função que o descreve quantitativamente é
a rugosidade, também denominada de largura de interface (“interface width”), que é
definida como [4]:
v
u
L
u1 X
W (L, t) = t
[h(i, t) − < h(i, t) >(t)]2
L i=1
(2.56)
Na equação 2.56, h(i, t) corresponde à altura da coluna num dado instante,
enquanto que < h(i, t) > (t) corresponde à sua média espacial, ou seja:
< h(i, t) > (t) =
L
1X
h(i, t)
L i=1
(2.57)
Se a taxa de deposição, i. e. o número de partı́culas que chegam em um sı́tio
por unidade de tempo, é constante, então podemos esperar que:
< h(i, t) > (t) ≈ t
(2.58)
O crescimento de uma superfı́cie como esta apresenta em geral duas etapas:
uma inicial, com forte variação de W (L, t), seguida de uma saturação do valor da
largura de interface, com (W (L, t) ≈ WSAT ). Um gráfico tı́pico da evolução temporal
da largura de interface possui portanto duas regiões distintas, separadas por um tempo
de saturação tx (“crossover time”). Isto sugere que a presença da saturação resulta da
existência de correlações no sistema, sendo um efeito inerente a sistemas finitos [4].
Para tempos menores que tx , W (L, t) aumenta com o tempo segundo uma lei
de potência de forma que:
W (L, t) ≈ tβ [t < tx ]
(2.59)
O expoente β, é chamado de expoente de crescimento e descreve a dependência
temporal do processo de formação da superfı́cie. Já para instantes maiores que tx ,
existe uma lei de potência relacionando a largura de interface e o tamanho do sistema
34
analisado, sendo então:
WSAT (L) ≈ Lα [t > tx ]
(2.60)
O expoente α, que caracteriza a rugosidade da interface saturada, é chamado
de expoente de rugosidade, e retrata como a rugosidade de um sistema é alterada com
a mudança da escala sob a qual o mesmo é analisado. O tempo de transição entre as
duas etapas também depende do tamanho do sistema, apresentando a seguinte relação
de escala:
tx (L) ≈ Lz
(2.61)
O expoente z, da relação acima, é chamado de expoente dinâmico. Os expoentes definidos não são independentes e, portanto, observa-se que no comportamento
da razão
W (L,t)
WSAT
com relação ao tempo, tem-se curvas que saturam ao mesmo tempo,
independentemente do tamanho do sistema. Analisando-se também o comportamento
t
,
tx
da largura de interface, como função de
mente. Estas observações sugerem que
tem-se curvas que saturam simultanea-
W (L,t)
WSAT
é uma função de
µ
W (L, t)
t
≈f
WSAT (L)
tx
t
tx
e portanto :
¶
(2.62)
A partir das equações 2.60, 2.61 e 2.62 obtém-se a relação de escala de FamilyVicsek [4], onde:
µ
t
W (L, t) ≈ L f
Lz
α
Fazendo a função de escala f
³
t
tx
´
¶
(2.63)
= f (u), verifica-se que, para pequenos
valores de u, a mesma cresce como uma lei de potência, de tal forma que f (u) ≈ uβ .
Para (t → ∞), a largura de interface satura, de modo que a função de escala é uma
constante para u À 1.
Para mostrar a relação entre os expoentes tratados, é importante observar
que sistemas de tamanhos diferentes possuem larguras de saturação diferentes, assim
35
Figura 2.13: Ilustração esquemática mostrando a variação da Largura de Interface com respeito
ao tempo para sistemas de diferentes tamanhos L, com L1 > L2 > L3 .
como tempos de saturação distintos (Figura 2.13). Analisando-se o comportamento
de log
³
W (L,t)
Lα
´
em função de log t, verifica-se que as curvas de saturação colapsarão
nas ordenadas para diferentes instantes tx (Figura 2.14). Finalmente, tratando o
comportamento de log
³
W (L,t)
Lα
´
em relação a log
³
t
tx
´
observa-se que as curvas de
saturação colapsarão, sugerindo a hipótese de escala (Figura 2.15).
Se o limite para o “ponto de crossover” (tx , W (tx )) for tomado pela esquerda,
tem-se, pela relação 2.59, que W (tx ) ≈ tβx . Por outro lado, se o limite, em relação ao
mesmo ponto, for tomado pela direita, tem-se, pela relação 2.60, que W (tx ) ≈ Lα .
Logo, pode-se concluir que, no ponto de “crossover”, tβx ≈ Lα e portanto:
Lz.β ≈ Lα ⇒ z =
α
β
(2.64)
Uma questão, que parece elementar, está relacionada à causa da saturação de
um sistema como este. O fato do tempo de saturação e a largura de interface saturada
crescerem com o tamanho do sistema, sugere que o fenômeno de saturação constitui
um efeito inerente a sistemas finitos. Uma vez que durante a deposição a superfı́cie
36
Figura 2.14: Representação logarı́tmica da divisão dos valores de rugosidade por Lα . Isto resulta
em que as curvas saturam para um mesmo valor de ordenada, porém em instantes distintos.
Figura 2.15: Comportamento logarı́tmico da rugosidade dividida por Lα e dos valores do tempo
por Lz . Neste caso saturam para um mesmo valor de ordenada e de abcissa.
37
obtida apresenta correlações, os diferentes sı́tios da superfı́cie formada por deposição
balı́stica não são completamente independentes, dependendo portanto da altura dos
sı́tios vizinhos.
Uma conseqüência é que as flutuações das alturas se propagam lateralmente,
uma vez que a próxima partı́cula depositada pode ter uma altura maior ou igual
à da sua vizinha, sendo portanto um processo de crescimento local. O crescimento
lateral, portanto, traz informações sobre como a altura de cada um dos vizinhos se
propaga globalmente. A distância tı́pica entre as alturas de cada sı́tio, em relação aos
vizinhos, ou seja, a distância caracterı́stica sobre a qual os sı́tios estão correlacionados,
é chamada de comprimento de correlação e é denotado por ξ|| . No começo do
crescimento os sı́tios são descorrelacionados, sendo que durante a deposição, ξ|| cresce
com o tempo inferindo-se que, para sistemas finitos, tal crescimento não é indefinido.
Quando ξ|| alcança o tamanho do sistema, a interface inteira se torna correlacionada,
resultando na saturação da largura de interface. Então tem-se que ξ|| ≈ L para
t À tx . De acordo com a relação 2.63, a saturação ocorre no instante tx dado por
2.61 e, portanto, relacionando o comprimento de correlação com o tempo, tem-se para
instantes muito menores que o “crossover time” que:
1
ξ|| ≈ t z
(2.65)
Uma equação contı́nua, chamada de equação KPZ [4], pode ser construı́da
envolvendo os expoentes caracterı́sticos deste modelo de deposição. Esta equação
constitui-se portanto uma classe de universalidade correspondente (ver Apêndice B).
38
Capı́tulo 3
Campo Elétrico Gerado por
Fronteiras Auto-Afins
Neste capı́tulo apresenta-se uma discussão sobre as propriedades de escala
de uma famı́lia de linhas e de superfı́cies equipotenciais numa região confinada entre
uma reta ou plano e um perfil ou seuperfı́cie com geometria fractal, região na qual
adota-se uma solução numérica para a equação de Laplace. Será feita a priori, uma
breve discussão a respeito das equações diferenciais parciais, bem como a respeito das
soluções analı́tica e numérica da equação de Laplace numa abordagem bi-dimensional.
Não obstante, a equação de Laplace também foi resolvida numericamente para sistemas
com 2D+1 dimensões, em uma região delimitada por uma fronteira fractal autoafim. Na primeira parte dos resultados, um estudo foi feito para o caso 1D+1,
investigando-se como o tamanho do sistema influencia as irregularidades das linhas
equipotenciais, ou seja, no valor das correspondentes dimensões fractais. Estas últimas
foram calculadas, para este caso, utilizando-se o método do Semivariograma. Na
segunda parte, foram investigadas as propriedades topológicas das superfı́cies equipotenciais (caso 2D+1), considerando-se superfı́cies base condutoras, geradas por uma
deposição balı́stica, bem como por uma metodologia Fractional Brownian Motion
(FBM). Neste caso foram utilizados os métodos Root Mean Square (RMS) e Semivariograma, possibilitando um estudo comparativo entre os dois métodos, bem como
inferindo os domı́nios de validade para os mesmos. É importante ressaltar que em todo
39
o trabalho, foram atribuı́das para as grandezas fı́sicas investigadas unidades efetivas.
3.1
Equações Diferenciais Parciais
Muitos problemas cientı́ficos requerem, para o seu tratamento adequado, uma
abordagem matemática no campo das equações diferenciais parciais. Um simples e
importante exemplo é um problema tı́pico de distribuição de temperatura em um
pedaço de metal, no qual certos pontos estão sujeitos a uma temperatura fixa. Assim,
num ponto genérico, a temperatura é uma função das coordenadas do ponto [29].
Fenômenos fı́sicos, que são descritos por funções de mais de uma variável
independente, podem resultar tratados por uma equação diferencial parcial, ou seja,
uma equação que envolve derivadas parciais. Sendo u uma função de duas variáveis
independentes, x e y, verifica-se que três derivadas parciais de segunda ordem podem
estar presentes na descrição matemática de um dado problema, ou seja:
∂ 2u ∂ 2u ∂2u
,
,
.
∂x2 ∂y 2 ∂x∂y
(3.1)
Dependendo do valor dos coeficientes dos termos das derivadas de segunda
ordem, as equações diferenciais parciais podem ser classificadas como elı́pticas, parabólicas ou hiperbólicas. A forma, as propriedades e a solução das equações diferenciais
dependem da natureza do problema fı́sico, bem como das condições de contorno a que
estão sujeitas. Em geral podemos escrever uma equação diferencial parcial de segunda
ordem como:
A
∂2u
∂ 2u
∂ 2u
∂u ∂u
+
B
+
C
+ D(x, y, u, , ) = 0
2
2
∂x
∂x∂y
∂y
∂x ∂y
(3.2)
A depender do valor de B 2 − 4AC, as equações serão caracterizadas como:
Elı́pticas → B 2 − 4AC < 0
(3.3)
P arabólicas → B 2 − 4AC = 0
(3.4)
40
Hiperbólicas → B 2 − 4AC > 0
(3.5)
Assim, para B = 0 e A = C = 1 a equação diferencial é sempre elı́ptica. Em
particular, se D = 0, a equação é chamada de equação de Laplace que, em coordenadas
cartesianas, assume a forma:
∂ 2u ∂2u
+
= 0 ⇒ ∆u = 0
∂x2 ∂y 2
(3.6)
Nesta seção busca-se a solução para o potencial elétrico no interior de uma
região desprovida de fontes, atribuindo-se as condições de Dirichlet (Apêndice C)
nos contornos que delimitam a região, inferior e superiormente. Lateralmente são
impostas condições periódicas. Soluções exatas são conhecidas em casos de contornos
com simetrias muito simples, como cı́rculos ou retângulos em duas dimensões, ou
esferas, cilindros ou paralelepı́pedos em três dimensões, com valores de u também
dados por funções muito simples sobre o contorno. Desta maneira, para contornos
com uma geometria não muito simples, é necessário recorrer a soluções numéricas.
Entender o comportamento do potencial e do campo elétrico próximo a superfı́cies condutoras é necessário para a adequada interpretação de vários fenômenos
presentes em diversas técnicas experimentais. A ionização do gás, pela ação do campo
no microscópio iônico de campo, bem como a evaporação de átomos e de ı́ons nos
processos de dessorção induzida, são bons exemplos [30].
3.2
Solução Analı́tica da Equação de Laplace
Os dispositivos elétricos funcionam baseados na ação de campos elétricos
produzidos por cargas elétricas e campos magnéticos produzidos por correntes elétricas
que, por sua vez, são constituı́das por cargas elétricas em movimento. Para entender
o funcionamento de dispositivos elétricos devemos estar aptos a quantitativamente
avaliar estes campos, nestes dispositivos e em torno deles. Isso é possı́vel se dominarmos técnicas que nos permitam determinar estas quantidades e, para facilitar
41
a visualização espacial, que possibilitem uma representação gráfica das grandezas
envolvidas, o que facilita a interpretação dos fenômenos fı́sicos envolvidos. Em outras
palavras, devemos ser capazes de produzir mapas de campos que descrevam o comportamento dos fenômenos elétricos. Estes mapas normalmente representam linhas de
fluxo, superfı́cies equipotenciais e distribuições de densidades de carga, possibilitando
a obtenção de informações a respeito da intensidade do campo, de diferenças de
potencial, de energia armazenada, de densidades de cargas, de densidades de correntes,
etc.
Para isso sabe-se que a Lei de Gauss, na forma integral, para um conjunto
discreto de cargas é dada por:
Z
E.η̂dS = 4π
X
qi
(3.7)
i
S
onde a soma é efetuada sobre todas as cargas no interior da superfı́cie fechada
S. Para uma distribuição contı́nua de carga, a lei de Gauss tem a forma:
Z
Z
E.η̂dS = 4π
S
ρ(r)dv
(3.8)
V
onde V é o volume limitado por S. Para obter-se a forma diferencial utiliza-se
o Teorema da Divergência onde, para um dado vetor A(r), definido em um volume V
limitado por uma superfı́cie fechada S, tem-se a seguinte relação:
Z
Z
A.η̂dS =
S
(5.A)dv
(3.9)
V
Portanto, aplicando-se o Teorema da Divervência à Lei de Gauss, ter-se-á:
Z
(5.E−4πρ)dv = 0
(3.10)
V
o que levará a:
5.E =4πρ
42
(3.11)
A igualdade 3.11 corresponde à Lei de Gauss na forma diferencial. Uma vez
que o campo eletrostático é conservativo:
5 × E =0
(3.12)
Como conseqüência direta de 3.12 tem-se que, neste caso, o campo elétrico
poderá ser escrito como o gradiente de uma função escalar, função esta denominada
de potencial elétrico, ou seja:
E=−5ϕ
(3.13)
Portanto, em regiões do espaço onde não se tem distribuição de carga, 5.E = 0
e, conseqüentemente, o potencial escalar satisfará a equação de Laplace, ou seja:
∆ϕ = 0
(3.14)
Para resolver analiticamente a equação 3.14, em coordenadas retangulares e
em duas dimensões, pode-se utilizar o método da separação de variáveis. Neste caso
assume-se que ϕ pode ser expresso como o produto de duas funções ξ e Ψ sendo
ξ = f (x) e Ψ = f (y). Então:
ϕ = ξ.Ψ.
(3.15)
Substituindo a equação 3.15 em 3.14 obter-se-á:
Ψ
d2 ξ
d2 Ψ
+
ξ
= 0.
dx2
dy 2
(3.16)
Dividindo-se 3.16 por ξ.Ψ, tem-se que:
1 d2 Ψ
1 d2 ξ
. 2+
= 0.
ξ dx
Ψ dy 2
(3.17)
Portanto, como a soma dos dois termos acima é uma constante, e cada variável
é independente uma da outra, cada termo deverá ser igual a uma constante. Desta
43
maneira pode-se escrever:
d2 ξ
= c21 .ξ
2
dx
(3.18)
d2 Ψ
= −c21 .Ψ = c22 .Ψ
dy 2
(3.19)
e
Logo, a solução geral da equação 3.14 será dada por:
φ = (K1 .ec1 x + K2 .e−c1 x )(K3 .ec2 y + K4 .e−c2 y )
(3.20)
onde K1 , K2 , K3 e K4 são constantes a determinar a partir das condições de
contorno especı́ficas do problema.
3.3
Solução Numérica em Duas Dimensões
O campo e o potencial elétrico são importantes quantidades na interpretação
de resultados de técnicas de imageamento de superfı́cies. Freqüentemente, a variação
local destas quantidades, em escala atômica, deve ser determinada com o fim de se
interpretar corretamente os dados experimentais. O efeito de uma única protuberância
em uma superfı́cie lisa tem sido convenientemente modelado a partir de caminhos
analı́ticos e numéricos. Porém, quando o perfil apresenta muitas irregularidades, o
problema se torna mais difı́cil, sendo portanto necessário recorrer a aproximações
numéricas para resolver a equação de Laplace. Esta equação apresenta uma solução
analı́tica para casos em que condições de contorno são aplicadas em sistemas fı́sicos de
alta simetria. Para resolvê-la numericamente, em uma região limitada por dois perfis
com potenciais constantes mas distintos e condições periódicas numa direção, usou-se
o método de Liebmann (Apêndice C). Tal método consiste em substituir derivadas
parciais por uma razão de diferenças. Torna-se necessário, para a implementação do
método, que o domı́nio, ou seja, a região onde se deseja determinar o potencial ponto
a ponto, seja divida em um “grid”. Valores iniciais são atribuı́dos arbitrariamente a
44
cada um dos pontos internos do “grid”, obedecendo à condição de que estes valores
estejam compreendidos entre os valores dos potenciais previamente atribuı́dos a dois
dos contornos que delimitam a região de integração. O potencial elétrico é então
calculado em cada nodo. Isto é feito considerando o laplaciano na forma:
∂ 2φ ∂ 2φ
+
=0
∂x2 ∂y 2
(3.21)
o que leva a:
∆φ =
φi+1,j − 2φi,j + φi−1,j φi,j+1 − 2φi,j + φi,j−1
+
(∆x)2
(∆y)2
(3.22)
Como ∆φ = 0, o potencial em cada ponto do “grid” ficará dado por:
φi,j =
sendo L =
∆x
∆y
φi+1,j + φi−1,j φi,j+1 + φi,j−1
+
2(L2 + 1)
2(L−2 + 1)
(3.23)
= 1, uma vez que neste trabalho optou-se por tratar o “grid”
com isometria na disposição dos nodos (“grid” quadrado).
3.4
Discussão e Resultados (1D+1)
Neste trabalho buscou-se considerar alguns aspectos adicionais a trabalhos
anteriores que tratam da influência de um condutor linear e rugoso, mantido a um
potencial constante e distinto do associado a um perfil retilı́neo distante [31, 32]. No
trabalho de Cajueiro et al. [31], o perfil utilizado foi um fractal auto-similar, a curva de
Koch. Posteriormente, Dias Filho et al. [32] consideraram perfis diferentes, resultado
de uma distribuição balı́stica e uma randômica com difusão, além de um perfil tipo
função de Weierstrass. A influência do tamanho do sistema não foi considerada nos
dois trabalhos imediatamente anteriores. Desta maneira, busca-se aqui discutir os
efeitos do tamanho do perfil na dimensão fractal de linhas equipotenciais compreendidas entre os dois perfis associados a potenciais fixos. Discute-se a rugosidade do perfil
e das linhas equipotenciais e sua dependência com o tamanho do perfil, com aplicação
a um caso de um perfil auto-afim [33].
45
Problemas deste tipo são motivados pela necessidade de análise do movimento
de partı́culas carregadas responsáveis pela formação de imagens em microscopias que
são resultado da trajetória de partı́culas sujeitas a ação de campos elétricos, como
é o caso das microscopias iônica de campo e de emissão de campo (FIM and FEM,
respectivamente). Efeitos resultantes das variações do campo elétrico local podem ser
também de interesse prático em outras áreas, como por exemplo para a determinação
das propriedades de dispositivos emissores [34], variação da função trabalho local e
comportamento de eletrodos [35, 36].
Para a determinação das propriedades do campo elétrico gerado por um condutor com geometria irregular, atribuiu-se ao mesmo um potencial elétrico constante.
Assim, no caso φ0 = 0. O perfil utilizado possui geometria invariante de escala,
correspondendo ao perfil mostrado na Figura 2.4. No nı́vel N1 , o perfil é constituı́do
por segmentos de reta que ligam sucessivamente os pontos representados por (x, y) =
(0, 0), (1, 1), (2, 0), (3, 1) e (4, 2).
Uma comparação das figuras em duas gerações
sucessivas, indica que todos os comprimentos horizontais foram ampliados por um
fator sh = 4, enquanto o fator de escala na direção vertical foi de sv = 2. Para
se evitar efeitos indesejáveis de fronteira, utilizou-se apenas metade dos pontos em
cada geração. Deste modo, para a metade da geração G8 , tem-se um perfil com 8192
pontos, conforme mostrado na Figura 3.1. Como o perfil possui caracterı́stica de autoafinidade, o expoente de rugosidade, α, considerando-se a construção já discutida,
possui valor 0, 5.
A outra condição de contorno está associada a uma linha reta condutora,
distante do perfil fractal, de modo que o potencial elétrico associado a este perfil, φ1 , foi
no caso tomado como igual a 100. A equação de Laplace foi resolvida no domı́nio entre
os dois condutores utilizando-se o método de Liebmann. O domı́nio foi convertido
em um “grid” bi-dimensional e o potencial elétrico foi calculado iterativamente em
cada ponto do mesmo, com condições de contorno periódicas impostas nas laterais
da região de integração. A partir da solução numérica da equação de Laplace, um
conjunto de linhas equipotenciais (Figura 3.2) cujas coordenadas são representadas
por yφi , i = 1, ..., M , foram determinadas por interpolação linear a partir dos valores
46
Figura 3.1: Porção de perfil auto-afim tratado como condutor, com 8192 pontos. A
Lei de Formação para este perfil é apresentada na seção 1.3.
dos potenciais do “grid”, φ.
As propriedades de rugosidade de cada perfil, yφi , seguem Leis de Escala
expressas pelos expoentes de rugosidade α, que foram calculados a partir do algoritmo
do semivariograma, conforme discutido no Capı́tulo 2. Como já se sabe, o mesmo é
baseado no cálculo da semivariança γφi (r) para cada perfil equipotencial, sendo esta
definida por:
X
1 n(r)
γφi (r) =
[yφi (x) − yφi (x + r)]2 .
2n(r) i=1
(3.24)
Na equação 3.24, yφi (x) indica a altura dos pontos que pertencem a um perfil
equipotencial e n(r) é o número de pares de pontos ao longo do perfil, separados
por uma distância r. Como a rugosidade do perfil, que mede essencialmente como
a diferença entre as alturas de dois pontos do perfil obedece a uma Lei de Escala, a
semivariança γ depende assintoticamente de r. Como visto no Capı́tulo 2 (Eq. 2.49),
pode-se concluir que:
47
Figura 3.2: Linhas equipotenciais φ(x) = φi , φi = 1, 7, 5n, n = 2, ..., 20, para o perfil
mostrado na Figura 3.1. φ = 100 para a fronteira superior é uma aproximação de
φ(y → ∞) → ∞.
γφi (r) ∼ r2αi
(3.25)
Para melhor entender os resultados para os perfis equipotenciais, divide-se a
discussão em duas partes. Primeiramente os resultados apresentados serão associados
a um perfil cujo número de pontos é mantido fixo. Na Figura 3.2 é mostrado um
conjunto de linhas equipotenciais. Cada linha equipotencial contém 8192 pontos,
correspondendo a um comprimento Lx = 8192 unidades de distância. Na direção
y tem-se 500 pontos, correspondendo a Ly = 500 unidades de comprimento. As 19
linhas mostradas na Figura 3.2 foram obtidas, para este caso especı́fico, considerandose o potencial elétrico do perfil base tal que φy1 (x) = 0 e o potencial da linha distante
φy2 (x) = 100. Os cálculos dos expoentes αi foram feitos para i = 1, ..., 100. Notase que, como as dimensões máximas, horizontal e vertical, são bem distintas (500 e
8192) a região onde a equação de Laplace foi resolvida é distorcida, no sentido de que
o comprimento do perfil é muito maior do que o máximo valor de y associado ao perfil
48
base (altura do perfil). Isto decorre da própria construção do perfil escolhido, onde a
distância máxima entre posições extremas na direção x cresce por um fator de escala
4, em cada passo da construção do perfil, enquanto na direção y está associado um
fator de escala 2. Isto leva ao fato de que a máxima altura do perfil cresce com o
1
comprimento Lx obedecendo à lei de potência L 2 .
A solução numérica da equação de Laplace, cuja representação gráfica é mostrada na Figura 3.2, mostra que a forma das linhas equipotenciais, partindo-se de uma
região próxima ao condutor base e em seguida daı́ se afastando, rapidamente perde
correlação com a forma do perfil.
Figura 3.3: (a)Semivariogramas para uma famı́lia de linhas equipotenciais calculadas
em uma região compreendida entre um perfil auto-afim tomado como base (Figura
3.1) e uma linha reta distante. Às linhas equipotenciais mais distantes do perfil
correspondem maiores inclinações do semivariograma no limite em que r → 0, o que
resulta em menores valores para a respectiva dimensão fractal. As curvas indicam,
de cima para baixo, as linhas equipotenciais cujos potenciais são φ = 1, 5 e 10n,
n = 1, 2, ..., 9. (b) Comportamento do desvio padrão relativo ao ajuste linear dos
semivariogramas mostrados em (a) como função do potencial elétrico, para a obtenção
dos correspondentes expoentes de rugosidade. O intervalo de escala utilizado para o
1
1
e Lrx ≤ 16
, onde Lrx representa a razão entre a escala
ajuste, foi tal que Lrx ≥ 800
correspondente e o tamanho do perfil.
Na Figura 3.3 mostram-se alguns semivariogramas tı́picos, em escala log-log,
para várias linhas equipotenciais com Lx = 8192. Vale observar que a maior escala
corresponde, aproximadamente, à metade do perfil inteiro. Cada ponto do semivari49
ograma representa o valor médio sobre todo o perfil, associado aos pares de pontos
separados por r. A forma de todos os semivariogramas, γ(r), em pequenas escalas,
corresponde à de pontos dispostos ao longo de uma linha reta em um gráfico loglog, com uma inclinação muito bem definida para um mesmo potencial, descrescendo
à medida que o potencial cresce. Porém, para escalas acima de 1000, a inclinação
decresce e o semivariograma γ(r) tem uma tendência à saturação, apresentando
aproximadamente a forma de um “plateau” horizontal. Isto decorre do fato de que,
para um perfil finito, uma escala maior que 1/4 do tamanho do perfil resulta em uma
estatı́stica “pobre” para o tratamento adequado do problema (ver Figura 3.4), ou seja,
existem poucos pares de pontos para a determinação da semivariância considerando-se
escalas maiores que 1/4 do tamanho do perfil. O mesmo efeito é observado para todos
os semivariogramas correspondentes às diferentes linhas equipotenciais.
Figura 3.4: Semivariogramas para dois perfis condutores com tamanhos Lx = 8192 e
Lx = 2048. Observa-se que para o perfil de menor tamanho, tem-se uma saturação
para menores comprimentos de escala em relação a um perfil de comprimento maior.
Este fato sugere que a saturação do semivariograma não depende de caracterı́sticas
topológicas do perfil.
É necessário então considerar a porção de cada semivariograma com comporta50
mento linear, de modo a se determinar o coeficiente angular de cada linha representada
na Figura 3.3, o que leva à determinação do expoente de rugosidade. Neste trabalho, o
coeficiente angular, que corresponde ao dobro do expoente de rugosidade α, foi obtido
pelo método dos mı́nimos quadrados. Para um valor fixo do comprimento do perfil,
Lx , o intervalo de escala utilizado para ajuste linear ficou compreendido entre
e
r
Lx
≤
1
.
16
r
Lx
≥
1
800
Naturalmente este intervalo está limitado à região anterior à saturação do
semivariograma. Os resultados indicam maiores valores de α para maiores valores de
φ. Portanto, o crescimento no valor de α implica no decréscimo do valor da dimensão
fractal Df . Assim, coerentemente com o que é possı́vel observar visualmente pelo
exame da Figura 3.2, nota-se que as linhas equipotenciais se tornam “menos rugosas”
quão mais distantes do perfil base.
Para avaliar a dependência da dimensão fractal com relação ao potencial elétrico
e/ou à distância média de cada linha equipotencial ao condutor base, consideramos
perfis de tamanhos Lx = 8192, Lx = 4096 e Lx = 2048. Na Figura 3.5, observa-se o
comportamento da Df − 1 como função do potencial elétrico φ, enquanto na Figura
3.6 Df − 1 é representada como função da distância média entre a linha equipotencial
e o perfil base, média esta definida por:
< d(φ) >=
L
1 X
[yφ (x) − yφ=0 (x)]
L x=0 i
(3.26)
Apesar da tendência geral de decaimento dos pontos nas Figuras 3.5 e 3.6,
nenhuma delas pode ser aproximada exatamente por exponenciais ou leis de potência.
Esta mesma caracterı́stica foi verificada para perfis randômicos [32] com os mesmos
expoentes de rugosidade α. Como, neste trabalho, considerou-se um perfil rugoso com
maior regularidade, i. e., com Lei de Formação definida, este resultado sugere uma
dependência não trivial entre a dimensão fractal Df e φ e/ou entre Df e < d(φ) >.
É interessante discutir o comportamento da distância média ao perfil base, de
cada linha equipotencial, com respeito ao potencial elétrico. Nota-se, como esperado,
que < d(φ) > cresce monotonicamente com o crescimento de φ. Contudo não foi
possı́vel encontrar uma relação trivial entre as duas quantidades. O exame detido
51
Figura 3.5: Comportamento da Df − 1 como uma função do potencial elétrico φ para
os três perfis. Cı́rculos, quadrados e triângulos indicam, respectivamente, perfis com
8192, 4096 e 2048 pontos.
Figura 3.6: Dependência da Df − 1 como função da distância média < d(φ) > com
respeito a φ. Cı́rculos, quadrados e triângulos indicam, respectivamente, perfis com
8192, 4096 e 2048 pontos.
52
da Figura 3.7 sugere a existência de duas regiões, uma para menores e outra para
maiores valores de φ, onde a dependência é aproximadamente linear. Para valores
intermediários esta simples dependência é perdida. Isto pode ser melhor analisado
mediante a avaliação do coeficiente de correlação e/ou do desvio padrão em grupos
de quatro pontos consecutivos do potencial mostrado na Figura 3.7. Para um dado
conjunto de pontos foram calculados os coeficientes angulares a partir do melhor
ajuste para quatro pontos sucessivos. Assim, o conjunto I corresponde ao ajuste dos
pontos 1, 2, 3 e 4, o conjunto II aos pontos 2, 3, 4 e 5 e assim sucessivamente. Por
fim, o conjunto VII corresponde ao ajuste linear dos pontos 7, 8, 9 e 10. Na Figura
3.8, mostra-se o desvio padrão referente ao ajuste linear para se obter os coeficientes
angulares, para cada conjunto de quatro pontos. Conforme se verifica, os resultados
são “mais pobres” para valores intermediários do potencial elétrico como discutido
antes. Para d e φ → ∞ o caráter de quase linearidade é recuperado.
Figura 3.7: Dependência da distância média < d(φ) > com respeito ao potencial
elétrico. Circunferências, quadrados e triângulos (os dois primeiros quase coincidentes)
indicam, respectivamente, perfis com 8192, 4096 e 2048 pontos.
Como segunda parte desta discussão, considerou-se também a dependência das
propriedades de rugosidade das linhas equipotenciais como uma função do comprimen53
Figura 3.8: Desvio padrão para um ajuste linear de um conjunto de quatro pontos
sucessivos. Os conjuntos são formados de modo que o conjunto I corresponde aos
pontos 1, 2, 3 e 4, o conjunto II corresponde aos pontos 2, 3, 4 e 5, e assim
sucessivamente.
to do perfil Lx . Este procedimento é útil no entendimento de como uma seqüência de
resultados, para amostras de comprimentos distintos mas todos finitos, pode indicar
o tipo de comportamento esperado no limite de um sistema de comprimento infinito.
Para isto será analisado, para cada linha equipotencial, como a dimensão fractal
associada a esta linha, varia com respeito a diferentes comprimentos Lx do perfil. Para
efetuar a comparação proposta, observou-se não ser necessário resolver seqüencialmente a equação de Laplace para sistemas de diferentes tamanhos Lx . Restrigiu-se assim
a integração da referida equação apenas para a região delimitada pelo perfil de maior
comprimento Lxmax , adotando-se, para comprimentos sucessivamente menores, os
valores do potencial segundo previamente calculado para um “grid” maior. Otimizase portanto o tempo computacional. Assim, é suficiente considerar, para cada linha
equipotencial, como a sua dimensão fractal depende de Lx , este considerado como
sendo de diferentes tamanhos. Adotando-se esta aproximação, não apenas o efeito
das condições periódicas é evitado, como também o problema de identificação das
54
linhas equipotenciais a serem estudadas, no caso de adotarem-se etapas distintas da
integração numérica, uma para cada tamanho de perfil. Como já visto, a determinação
da dimensão fractal requer uma prévia determinação do expoente de rugosidade. É
necessário portanto discutir um aspecto importante que está associado à medida
do expoente de rugosidade para perfis de diferentes comprimentos, qual seja uma
adequada adoção das escalas inicial e final para o ajuste em busca da linearidade
conforme já descrito. Com o crescimento de Lx , é observado que o intervalo de
escala r utilizado nos semivariogramas também pode ser maior, de modo que os
valores dos expoentes de rugosidade dependerão dos maiores comprimentos de escala
a serem incluı́dos no processo de linearização. Conforme já mencionado, a escala
máxima utilizada foi rmax =
Lx
.
16
Portanto utilizaram-se, no estudo da presente
seção, dois caminhos distintos para a interpretação dos resultados. Dois resultados
aparentemente contraditórios foram encontrados: i) a dimensão fractal das linhas
equipotenciais de um perfil com comprimento fixo decresce quando sucessivamente
consideramos linhas mais afastadas do perfil base; ii) a dimensão fractal de uma
determinada linha equipotencial (valor de φ mantido fixo) cresce com o aumento do
tamanho do perfil (Figura 3.5). Este segundo resultado, sugere que, tendo-se um perfil
infinito, todas as linhas equipotenciais serão caracterizadas pelos mesmos valores dos
expoentes de rugosidade apresentando portanto as mesmas dimensões fractais. Vale
observar que a interpretação do segundo resultado se deve a propriedades de escala
das linhas equipotenciais, com as conclusões para o sistema hipoteticamente infinito
não apresentando contradição com o caso finito.
3.5
Discussão e Resultados (2D+1)
Nesta seção, é feita uma discussão a respeito das propriedades fractais de
superfı́cies equipotenciais.
Para isto foi imprescindı́vel uma extensão da solução
numérica da equação de Laplace para o caso de 2D + 1 dimensões. Os resultados
obtidos tiveram o seu tratamento dividido em duas partes: primeiramente, discutese algumas propriedades das superfı́cies equipotenciais calculadas para condutores
55
base obtidos por deposição balı́stica e por uma metodologia FBM. Como segunda
parte, uma discussão complementar dos resultados é feita com respeito ao domı́nio de
validade dos métodos utilizados para quantificar as irregularidades que caracterizam
as superfı́cies equipotenciais.
A priori, o algoritmo para a solução numérica da equação de Laplace utilizado
no estudo de perfis fractais (Caso 1D+1), foi extendido para o caso de 2D + 1
dimensões (Apêndice C). Isto possibilitou o cálculo das superfı́cies equipotenciais
numa região limitada por uma superfı́cie resultado de uma deposição balı́stica DB ou
de uma metodologia FBM, e uma superfı́cie plana distante. A diferença de potencial
elétrico entre as duas superfı́cies foi assumida constante, de modo que para o condutor
base atribuiu-se o potencial elétrico φ0 = 100 e para o condutor distante φ1 = 0,
diferentemente portanto, das atribuições adotadas no estudo feito para o caso 1D+1.
Através do método de Liebmann, com isometria na disposição dos pontos do “grid”,
extendeu-se a expressão para o cálculo do potencial elétrico para pontos do espaço
(Apêndice C). Uma vez determinado o potencial elétrico nos pontos do “grid”, um
algoritmo foi construı́do no sentido de se determinar o potencial elétrico associado a
qualquer ponto do espaço, utilizando-se para isso uma interpolação linear. Definindose a variação do potencial com respeito à direção definida pelo versor k̂ (versor que
representa a direção do eixo que define a altura de cada ponto do espaço.) e denotandose tal por Gφ , então tem-se que:
Gφ =
φ(i, j, k + 1) − φ(i, j, k)
∆k
(3.27)
com ∆k = k + 1 − k = 1, uma vez que o espaçamento dos pontos vizinhos do
“grid” foi considerado unitário ( Um estudo envolvendo o processo de relaxação para a
convergência obtida na solução numérica da equação de Laplace pode ser encontrado
em [31] ) . As coordenadas de uma superfı́cie equipotencial qualquer são especificadas
por (i, j, k + dz), com 1 > dz > 0, de modo que o potencial elétrico que se busca será
dado por:
φ(i, j, k + dz) = φ(i, j, k) + Gφ .dz
56
(3.28)
Uma vez que dz é perfeitamente determinado pela relação 3.28, as coordenadas
da superfı́cie equipotencial que se busca serão completamente especificadas.
Figura 3.9: Superfı́cies equipotenciais calculadas para o caso de uma superfı́cie base irregular
gerada por uma metodologia FBM, com tamanho 200x200.
No caso dos resultados obtidos para uma superfı́cie condutora obtida por DB,
duas situações foram consideradas. Uma vez que o comprimento de correlação lateral
é não nulo durante o processo de deposição, após a saturação da rugosidade verificase a presença de vazios no volume (“bulk”) do material. Este fato acarreta um
aumento do tempo computacional, uma vez que o processo iterativo para solução
da equação de Laplace será considerado para tais pontos que pertencem aos vazios.
Para a “construção” da superfı́cie uma aproximação foi necessário ser feita. Uma
vez que na deposição balı́stica podem surgir vazios em posições abaixo de pontos
que pertencem ao material condutor, nestes vazios, teria-se pontos cujo potencial
elétrico seria determinado a partir do processo de integração numérica. Isto traria
conseqüências que resultariam na existência de uma superfı́cie externa e de uma
superfı́cie interna, relativamente à porção de espaço resultante da deposição. Assim
optou-se por ignorar as porções correspondentes aos vazios, associando o potencial
dos seus pontos interiores com o potencial da fronteira. Portanto, considerando-se o
potencial elétrico da superfı́cie como φSup (x, y, zSup ), foi feita a aproximação:
57
φ(x, y, z < zSup ) = φSup (x, y, zSup )
(3.29)
Houve também uma preocupação durante o processo de deposição balı́stica
que está associada à saturação da rugosidade da superfı́cie condutora. Para isso,
foi considerado o número de partı́culas a serem depositadas (ver Figura 3.10 ) para
a configuração da superfı́cie saturada. Este número naturalmente aumenta com o
tamanho do sistema, garantindo portanto uma superfı́cie com rugosidade definida.
Isto ganha uma importância tecnológica, pois em se tratando de um catalisador
percebe-se no processo de sı́ntese que a partir de uma determinada quantidade de
partı́culas depositadas, ocorre a saturação dos valores de dimensões fractais associadas
ao material no suporte. Nesse ponto, a irregularidade da superfı́cie é máxima e,
portanto, um aumento da quantidade de material depositada no suporte não significa
melhoria da qualidade do catalisador [37] .
Figura 3.10: O gráfico acima mostra o comportamento da rugosidade com relação ao tempo para
superfı́cies obtidas por deposição balı́stica com tamanhos 100x100, 200x200 e 300x300. Observa-se
um crescimento preliminar da rugosidade, para uma posterior saturação.
Para quantificação das irregularidades das superfı́cies equipotenciais utilizaramse os Métodos RMS e Semivariograma discutidos no Capı́tulo 2. Nas Figuras 3.11 e
58
3.12, verifica-se o comportamento do desvio quadrático médio (RMS) como função
do comprimento de escala utilizado, em escala logarı́tmica. O cálculo foi feito para
superfı́cies equipotenciais φi = 99 − 3n, n = 0, ..., 6, em que os condutores base
foram obtidos por DB, bem como por uma metodologia FBM. Foram considerados,
sistemas cujos tamanhos são 200x200 e 300x300. Os expoentes de rugosidade foram
determinados tomando como escalas inicial e final 2 e 20 respectivamente, para o
caso de um sistema 200x200 e 2 e 30 para o caso de um sistema 300x300. Para
um sistema 200x200x100, bem como 300x300x100 de volume, verifica-se claramente
que, para pequenas escalas, a diferença entre as rugosidades das correspondentes
superfı́cies equipotenciais é maior, observando-se uma diminuição dessa diferença
para maiores comprimentos de escala. Isto pode ser justificado pelo fato de que
superfı́cies equipotenciais mais próximas do condutor base são caracterizadas por
menores expoentes de rugosidade apresentando, conseqüentemente, maiores valores de
dimenão fractal, resultado esse em concordância com o estudo feito para perfis fractais.
A Figura 3.13, mostra o comportamento do desvio padrão obtido para o ajuste linear
dos gráficos 3.11 e 3.12 como função do potencial elétrico, para a determinação dos
expoentes de rugosidade de acordo com os intervalos de escala já especificados.
Outro importante aspecto nessa discussão está associado ao comportamento
da dimensão fractal com respeito ao potencial elétrico para os sistemas descritos
anteriormente (ver Figura 3.14). Para um sistema 200x200 as superfı́cies condutoras,
obtidas por DB e pela metodologia FBM, possuem rugosidades 5, 37 e 8, 64 respectivamente. Para um sistema 300x300 as correspondentes rugosidades obtidas foram
6, 49 e 5, 55.
A partir da Figura 3.14, observa-se que a geometria do condutor base afeta
claramente o grau de correlação geométrica entre o mesmo e as superfı́cies equipotenciais na sua imediata vizinhança. Comparando-se o comportamento da dimensão fractal
como função do potencial elétrico, para sistemas diferentes e de mesmo tamanho
verifica-se que, para o caso do condutor base obtido por DB, o decréscimo da dimensão
fractal com o decréscimo do potencial elétrico é menos acentuado, ou seja, a geometria
do condutor base é refletida em superfı́cies equipotenciais mais distantes do condutor
59
Figura 3.11: Comportamento do desvio quadrático médio com o correspondente comprimento
de escala no método RMS. As superfı́cies base consideradas neste caso foram obtidas por uma
metodologia FBM e por deposição balı́stica com tamanhos 200x200. Maiores expoentes de rugosidade
foram obtidos para menores valores do potencial elétrico. As curvas indicam, de baixo para cima, as
superfı́cies equipotenciais cujos potenciais são 99 − 3n, n = 6, 5, ..., 1, 0. Vale salientar que a unidade
do RMS é dada em unidades arbitrárias, assim como a escala.
Figura 3.12: Comportamento do desvio quadrático médio com o correspondente comprimento
de escala no método RMS. As superfı́cies base consideradas neste caso foram obtidas por uma
metodologia FBM e por deposição balı́stica com tamanhos 300x300. Maiores expoentes de rugosidade
foram obtidos para menores valores do potencial elétrico. As curvas indicam, de baixo para cima, as
superfı́cies equipotenciais cujos potenciais são 99 − 3n, n = 6, 5, ..., 1, 0. Vale salientar que a unidade
do RMS é dada em unidades arbitrárias, assim como a escala.
60
Figura 3.13: (a) Comportamento do desvio padrão como função do potencial elétrico no ajuste
linear das curvas mostradas na Figura 3.11 para o intervalo de escala 2 ≤ ω ≤ 20. (b) Comportamento
do desvio padrão como função do potencial elétrico no ajuste linear das curvas mostradas na Figura
3.12 para o intervalo de escala 2 ≤ ω ≤ 30.
base, de maneira mais significativa.
Este aspecto ganha importância quando se
investiga o efeito de várias protuberâncias que pertencem a uma superfı́cie mais ou
menos suave, na distribuição do campo elétrico no espaço. Os resultados indicam
que superfı́cies que são menos suaves, no caso as formadas por DB, contribuem mais
significativamente na intensidade do campo elétrico médio associado a superfı́cies
equipotenciais mais distantes. Este aspecto voltará a ser discutido com mais detalhes
no Capı́tulo 4 no tratamento da dinâmica de partı́culas carregadas, não interagentes,
neste campo irregular.
Considerando-se o método RMS observa-se claramente que, para valores de
potenciais próximos a 0, os valores das dimensões fractais das superfı́cies equipotenciais
correspondentes tendem a 2, 1. Isto é justificado pela sensibilidade do método em se
analisar superfı́cies rugosas. Uma vez que o desvio quadrático médio para um dado
comprimento de escala e para uma amostra analisada da superfı́cie é calculado, todos
os pontos para o intervalo de amostra estudado são considerados.
Comportamentos semelhantes são verificados do ponto de vista qualitativo,
utilizando-se o método semivariograma para a determinação das dimensões fractais
das superfı́cies equipotenciais. Na Figura 3.15, o comportamento da dimensão fractal
61
Figura 3.14: Os gráficos acima mostram o comportamento da dimensão fractal (Df) como função
do potencial elétrico φ, para superfı́cies obtidas por deposição balı́stica e por uma metodologia FBM
com tamanhos 200x200 e 300x300. Utilizou-se o método RMS para a determinação dos expoentes
de rugosidade.
das superfı́cies equipotenciais com o potencial elétrico é mostrado utilizando-se o
método semivariograma.
As tendências de decaimento da dimensão fractal com o decréscimo do potencial
(para os valores adotados de potencial na fronteira) apresentadas nas Figuras 3.14 e
3.15 para o caso dos condutores base estudados confirmam, como exposto na seção
anterior, uma não linearidade que não apresenta caracterı́sticas de funções exponenciais ou leis de potência.
Outro fato a observar é que, quantitativamente, os dois métodos apresentam
tendências diferentes, ou seja, as dimensões fractais apresentam valores distintos para
superfı́cies equipotencias mais próximas do condutor cuja geometria é plana. No
caso do método RMS verifica-se que o expoente de rugosidade tende a 0,9 (dimensão
fractal tende a 2,1) enquanto no método semivariograma essa tendência vai para 1,0
(dimensão fractal tende a 2,0) justificando-se então uma maior precisão do método
RMS, comparativamente ao do Semivariograma. (ver Figura 3.16).
Contudo, uma situação foi verificada durante a realização dos cálculos para
determinação dos expoentes de rugosidade a partir dos métodos RMS e Semivariograma. Os mesmos apresentaram uma restrição para o cálculo das dimensões fractais no
limite em que as superfı́cies equipotenciais apresentam geometrias próximas à de um
62
Figura 3.15: Os gráficos acima mostram o comportamento da dimensão fractal (Df) como função
do potencial elétrico φ, para superfı́cies obtidas por deposição balı́stica e por uma metodologia FBM
com tamanhos 200x200 e 300x300. Utilizou-se o método semivariograma para a determinação dos
expoentes de rugosidade.
Figura 3.16: Os gráficos acima mostram o comportamento dos expoentes de rugosidade com relação
ao potencial elétrico φ, para superfı́cies obtidas por deposição balı́stica e por uma metodologia
FBM com tamanhos 200 x 200. Verifica-se que α−→0, 9 no método RMS e α−→1, 0 no método
Semivariograma, inferindo-se uma menor precisão deste último. As barras de erro mostradas, estão
associadas ao ajuste linear para determinação dos expoentes de rugosidade.
63
plano, ou seja, apresentam expoentes de rugosidades próximos a 1, 0. Encontrou-se
uma dificuldade em se obter mais pontos para a caracterização do comportamento
da dimensão fractal com relação ao potencial elétrico para sistemas cujos tamanhos
nas direções x e y são da ordem ou menores que a altura da caixa que contém a
região de integração numérica e os contornos associados às condições de fronteira. Na
Figura 3.17 é mostrado o comportamento do expoente de rugosidade com relação à
rugosidade das superfı́cies equipotenciais analisadas, para o caso de um condutor base
gerado por deposição balı́stica. Os sistemas analisados foram tais que apresentavam
volumes 100x100x100, 200x200x100 e 300x300x100.
Figura 3.17: (a) Gráfico mostrando o comportamento do expoente de rugosidade, calculado
pelo método do Semivariograma, com relação à rugosidade para uma superfı́cie base gerada por
deposição balı́stica. (b) Ampliação da região do gráfico (a), em que a rugosidade tende a 0 e
o expoente de rugosidade tende a 1,0. Os pontos indicados por P L, correspondem à superfı́cie
equipotencial investigada que ainda apresenta propriedades de auto-afinidade, sendo portanto o
método semivariograma aplicável. As barras de erro estão associadas ao erro no cálculo da inclinação
do semivariograma para a determinação dos expoentes de rugosidade.
O ponto P L, indicado na Figura 3.17 representa a superfı́cie equipotencial
limite em que o método ainda teve validade. O critério utilizado para esta avaliação,
está associado à determinação dos expoentes de rugosidade para estas superfı́cies
equipotenciais, que apresentaram valores não convenientes com o esperado. Este
critério será discutido mais adiante, ainda nesta seção. Os valores de rugosidade
associados a tais superfı́cies, para os sistemas 100x100x100, 200x200x100 e 300x300x100, foram respectivamente 0, 083, 0, 079 e 0, 257. Já os expoentes de rugosidade cor64
respondentes foram 0, 9849, 0, 9780 e 0, 9750. Isto sugere que existe um expoente de
rugosidade médio limite, em torno do qual o método semivariograma não é aplicável.
Isto quer dizer que o correspondente método determina valores de dimensões fractais
que não são compatı́veis com os valores esperados para superfı́cies pouco rugosas.
Este efeito de limite de validade também foi verificado para o método RMS,
utilizando-se os mesmos sistemas 100x100x100, 200x200x100 e 300x300x100. Contudo
a maior precisão deste método, reflete expoentes de rugosidade limites diferentes. A
Figura 3.18 mostra o comportamento da rugosidade com respeito ao expoente de
rugosidade, utilizando-se o método RMS. Para este caso, os expoentes de rugosidade
limite correspondentes aos sistemas de volumes 100x100x100, 200x200x100 e 300x300x100 foram respectivamente 0, 8840, 0, 8835 e 0, 8905, com rugosidades já mencionadas no parágrafo anterior. Logo, a hipótese de que existe um expoente de rugosidade
médio limite a partir do qual o as superfı́cies equipotenciais perdem as caracterı́sticas
de auto-afinidade é confirmada também para o método RMS.
Figura 3.18: (a) Gráfico mostrando o comportamento do expoente de rugosidade calculado pelo
método RMS, com relação à rugosidade para uma superfı́cie base gerada por deposição balı́stica.
(b)Ampliação da região do gráfico (a), em que a rugosidade tende a 0 e o expoente de rugosidade
tende a 1,0. Os pontos indicados por P L, correspondem à superfı́cie equipotencial investigada que
ainda apresenta propriedades de auto-afinidade, sendo portanto o método RMS aplicável. As barras
de erro estão associadas ao erro no cálculo da inclinação do RMS para a determinação dos expoentes
de rugosidade.
A justificativa para a existência de domı́nios de validade nos métodos RMS
e Semivariograma está associada ao fato de que, para superfı́cies com geometrias
65
próximas à geometria euclidiana, as variações tanto do desvio quadrático médio quanto
da semivariância, são muito pequenas para grandes variações no comprimento de
escala a partir do qual a superfı́cie é caracterizada. Deste modo, a perda de caracterı́sticas fractais fica evidente, uma vez que as leis de potência que relacionam
o desvio quadrático médio e a semivariância com o expoente de rugosidade não
são mais consideradas. A Figura 3.19 retrata em (a), o comportamento do desvio
quadrático médio, com valores normalizados, como função da escala utilizada para
três superfı́cies equipotenciais calculadas, considerando-se o condutor base gerado por
DB com tamanho 100 x 100. Observa-se uma incompatibilidade na determinação do
valor do expoente de rugosidade para φ = 5, utilizando-se o método RMS, uma vez
que para este potencial elétrico tem-se um expoente de rugosidade menor em relação
a φ = 25. O comportamento de lei de potência é comprometido neste caso, uma vez
que o erro apresentado no ajuste linear é maior para φ = 5. Em (b) é mostrado o
comportamento do desvio quadrático médio como função da escala para as superfı́cies
equipotenciais consideradas em (a). Observa-se claramente pequenas variações do
desvio quadrático médio para grandes variações de escala, justificando a inabilidade
de utilização do método RMS para estes casos.
Por fim, nesta subseção, é apresentado o comportamento da distância média
das superfı́cies equipotenciais ao condutor base como função do potencial elétrico. A
distância média de uma determinada superfı́cie equipotencial ao condutor base, para
um sistema de tamanho L x L, é definida pela expressão:
L2
1 X
hd(φ)i = 2
[zφ (x, y) − zφ=100 (x, y)]
L n=1 i
(3.30)
A partir da Figura 3.20, observa-se um comportamento da distância média
com respeito ao potencial elétrico em concordância com os resultados obtidos no caso
1D+1 para potenciais intermediários.
Porém, uma situação complementar foi verificada. Observa-se que, no limite
em que a distância média das superfı́cies equipotenciais ao condutor base tende a zero,
seja por DB ou por uma metodologia FBM, um comportamento não linear é verificado
mais acentuadamente. Esse resultado é complementar aos resultados expostos para
66
Figura 3.19: O gráfico (a) acima, mostra o comportamento do desvio quadrático médio como
função da escala, em escala logarı́tmica para as superfı́cies equipotenciais φ = 70, φ = 25 e φ = 5. Os
expoentes de rugosidade bem como os erros associados ao cálculo também são mostrados. O gráfico
reduzido (b) mostra o comportamento do desvio quadrático médio como função da escala retratando
pequenas variações para φ = 5
Figura 3.20: Os gráficos acima mostram o comportamento da distância média de cada superfı́cie
equipotencial ao condutor base, como função do potencial elétrico. Estes cálculos foram feitos para
sistemas 200x200 e 300x300 para condutores base formados por DB e pela metodologia FBM.
67
o caso 1D+1, uma vez que no atual estudo, algumas superfı́cies equipotenciais com
distâncias médias ainda menores em relação ao condutor base foram consideradas.
Este aspecto ganha importância, no sentido de se investigar como, na média, a
geometria do condutor base afeta propriedades de rugosidade das superfı́cies equipotenciais próximas ao condutor base. A Figura 3.21, mostra este comportamento
não trivial, no limite em que < d(φ) >→0.
Figura 3.21: Comportamento da distância média normalizada de cada superfı́cie equipotencial ao
condutor base, como função do potencial elétrico para uma situação em que < d(φ) >→0. Foram
consideradas superfı́cies equipotenciais com potenciais 99, 1 ≤ φ ≤ 99, 9.
Foram considerados neste último, superfı́cies equipotenciais elétricas localizadas entre φ=99 e φ=100. Esta transição da não linearidade para situações onde o
comportamento de < d(φ) > com φ é mais aproximadamente linear, sugere que as
variações das rugosidades das superfı́cies equipotenciais próximas ao condutor base
com o potencial elétrico são mais acentuadas para pequenas variações do potencial
elétrico. Este fato pode ser verificado a partir da Figura 3.22. Para distâncias médias
intermediárias, a variação da rugosidade com o potencial elétrico é mais suave. Estes
resultados indicam portanto, que a geometria associada às superfı́cies equipotenciais
próximas ao condutor base, guardam informações locais mais significantes do condutor
base quando a essas regiões estão associados campos elétricos mais intensos.
68
Figura 3.22: Comportamento da rugosidade como função do potencial elétrico considerando-se
como condutor base o gerado por DB - (fig.esquerda) e por uma metodologia FBM - (fig. direita).
Constata-se variações acentuadas na rugosidade para pequenas variações do potencial elétrico,
considerando-se as superfı́cies equipotenciais próximas ao condutor base (99, 1 ≤ φ ≤ 99, 9). Para
potenciais elétricos intermediários a variação da rugosidade com a variação do potencial elétrico é
mais suave.
69
Capı́tulo 4
Dinâmica de Partı́culas Carregadas
em um Campo Elétrico Irregular
Neste Capı́tulo é estudada a dinâmica de partı́culas carregadas, não interagentes entre si, em um campo elétrico produzido por um condutor rugoso e com geometria
fractal. Será feita preliminarmente uma discussão a respeito dos problemas clássicos
envolvendo a dinâmica molecular, apresentando-se o método Preditor-Corretor. Posteriormente, na seção Discussão e Resultados, será tratado um algoritmo para a
determinação do campo elétrico em um ponto qualquer do espaço bem como um
estudo do comportamento médio deste campo para sistemas condutores de tamanhos
100x100. Situações em que os condutores foram gerados por DB e por uma metodologia
FBM foram analisadas. Por fim, será discutido como grandezas fı́sicas associadas à
dinâmica das partı́culas carregadas podem trazer informações a respeito da distribuição
espacial do campo elétrostático, bem como da geometria do condutor base.
4.1
Equações do Movimento
Com a evolução da simulação computacional ocorrida nos últimos vinte anos,
a Dinâmica Molecular tornou-se um método poderoso para o estudo da matéria condensada. Áreas de aplicação como a modelagem molecular, modelos de spin, lı́quidos
e a dinâmica molecular quântica podem ser citadas . Uma adequada formulação do
70
modelo a ser utilizado para a descrição da matéria é de fundamental importância no
tratamento computacional do sistema de interesse.
A dinâmica de um sistema clássico, não relativista, de N partı́culas interagentes
entre si, é descrita pelo formalismo Hamiltoniano por 6N equações do movimento, de
modo que as taxas de variação com respeito ao tempo, da posição e do momento
generalizado, são dadas respectivamente por:
q˙j =
∂H
∂pj
(4.1)
e
p˙j = −
∂H
∂qj
(4.2)
com j = 1, 2, ..., N . O termo H é o Hamiltoniano do sistema sendo função
das posições e dos momentos generalizados, denotados por q e p respectivamente, ou
seja, H = H(q, p, t). As equações do movimento formam um sistema de equações
diferenciais de primeira ordem, geralmente acopladas, que são resolvidas a priori
conjuntamente. Dada a impossibilidade de obtenção analı́tica destas soluções, a única
forma de se obter a evolução do sistema em estudo se dá através da implementação
de métodos numéricos e portanto aproximados.
É importante recordar que habitualmente um problema de dinâmica de partı́culas é tratado como um problema de Cauchy, ou seja, as condições iniciais que
definem a posição e a velocidade, ~q(t0 ) = ~q0 e p~(t0 ) = p~0 , são conhecidas. Assim,
existirá uma solução para as equações do movimento e a mesma será única.
A um sistema Hamiltoniano 6N dimensional associa-se um espaço de fase,
Γ, também 6N dimensional, em que as coordenadas são 3N posições generalizadas,
qj , e os 3N momentos generalizados, pj . O estado deste sistema será perfeitamente
determinado por suas 6N coordenadas no espaço de fase, de modo que a evolução do
sistema tratado dará lugar a uma trajetória neste espaço.
Neste trabalho o sistema estudado compõe os movimentos de partı́culas carregadas não interagentes entre si, sujeitas a um campo elétrico irregular. Para este
71
sistema, o hamiltoniano, independente do tempo, será dado por:
H(q, p) =
N
X
p2j
j=1
2mj
+
N
X
U (qj )
(4.3)
j=1
onde U (qj ) corresponde à energia de interação da partı́cula j com o campo
elétrostático externo. A equação do movimento para uma partı́cuja j com carga q e
massa m, interagindo com este campo elétrico irregular é dada por:
−
→
d2 →
rj
−q∇φ(−
rj )
=
2
dt
m
(4.4)
→
, sendo φ(−
rj ) a função potencial elétrico.
4.1.1
Método Preditor-Corretor
O fundamento para a simulação mediante a técnica de Dinâmica Molecular
(DM), termo como é conhecido na literatura, é o conhecimento das equações de
movimento para as partı́culas constituintes do sistema considerado. Por sua vez,
o algoritmo de um programa para a DM consiste em uma metodologia para a solução
numérica destas equações, fornecendo uma trajetória, i. e., coordenadas e momentos
conjugados em função do tempo, para as partı́culas do sistema em estudo. A solução
das equações, uma vez que numérica, implica em tratar o tempo descontinuamente,
com os intervalos entre um instante e o seguinte sendo denominado de passo de tempo
do processo integrativo. A escolha do passo depende do problema especı́fico que está
sendo tratado, podendo assim ser adaptado segundo as caracterı́sticas do problema
especı́fico.
Assim, a partir da trajetória, propriedades de equilı́brio e grandezas
dinâmicas podem ser calculadas.
Os algoritmos de integração temporal se baseiam em métodos de diferenças
finitas. Conhecendo as posições e algumas de suas derivadas em um instante t, o
esquema de integração proporciona as mesmas quantidades em um tempo posterior
t+ ∆t. Iterando este procedimento, se pode seguir a evolução destes sistemas ao longo
do tempo.
Aos esquemas de integração estão associados dois tipos de erro, que são i72
nevitáveis e competem entre si. No primeiro caso temos os erros de truncamento,
que são devidos à implementação de um algoritmo numérico para se resolver um
problema analı́tico, i. e., trata-se uma variável essencialmente contı́nua - o tempo
- como uma variável discreta. Em geral decorrem do desenvolvimento em séries de
Taylor, truncadas a partir de uma certa ordem. São portanto erros intrı́nsecos ao
algoritmo e não dependem da implementação computacional adotada. Outro tipo
de erro é o de arredondamento, que está associado ao fato de se trabalhar com uma
aritmética de precisão finita, caracterı́stica do cálculo computacional. Neste caso temse, por exemplo, o número finito de dı́gitos usados na representação de um número
real irracional.
Em geral pode-se dizer, heuristicamente, que os erros de truncamento podem
ser reduzidos pela escolha de um passo de tempo menor, com o que se teria uma maior
aproximação de uma variável contı́nua. Contudo, adotando-se este procedimento,
os erros de arredondamento serão aumentados consideravelmente, uma vez que a
quantidade de iterações feitas, para um mesmo intervalo total de tempo, crescem.
Por outro lado, um passo de tempo excessivamente grande diminuirá os erros de
arredondamento, aumentando-se portanto os erros de truncamento. A escolha do
passo levará portanto em conta o estado do sistema, uma vez que utilizando-se
de intervalos de tempo demasiadamente grandes resultará em uma descrição fı́sica
incompatı́vel com uma correta evolução do sistema. Da competição entre os dois tipos
de erro resulta a existência de um intervalo de tempo ótimo, em que a combinação
dos erros é minimizada de modo a se obter melhores resultados. Porém não existem
métodos analı́ticos que permitam obter um valor “ideal” para o passo.
Um tipo de algoritmo de integração para as equações de movimento é o denominado preditor − corretor [38]. Este método retém termos além dos de segunda
ordem na expansão de Taylor. Assim sendo, diferentemente do caso do algoritmo de
Verlet (ver Apêndice D), neste caso não se faz necessário a adoção de intervalos de
tempo tão pequenos. Além disso, este método tem a vantagem de permitir que se use
um intervalo de tempo variável.
A estrutura do algoritmo do preditor−corretor compreende a determinação da
73
posição, e das suas derivadas, em um instante tempo t + ∆t, a partir do conhecimento
das posições e das suas n derivadas num instante de tempo t. Este conjunto de
valores, para a primitiva e derivadas, constitui o que se denomina de valores preditos.
Posteriormente, com o conhecimento da posição predita, calcula-se a intensidade
da força atuante em cada partı́cula e portanto as acelerações, que correspondem
às derivadas de segunda ordem. A seguir, comparam-se as acelerações predita e a
determinada a partir da força, determinando-se assim o valor de uma correção para
este termo. Com isto corrigem-se as posições e suas derivadas, usando as correções
calculadas. Os novos termos para a primitiva e suas derivadas constituem os termos
corrigidos.
No presente caso utilizou-se a expansão, até a quinta ordem no intervalo de
tempo. Desta maneira, a expansão em série de Taylor da posição r~j de uma j-ésima
partı́cula pode ser expressa da seguinte forma:
5
→
X
dn −
rj (t) n
−
→
→
−
rj (t + ∆t) = rj (t) +
∆t
dtn
n=1
(4.5)
Por comodidade para a notação definiremos os termos da equação 4.5 como
sendo:
→
d−
rj
∆t
dt
(4.6)
βj (t) =
→
d2 −
rj
(∆t)2
dt2
(4.7)
γj (t) =
→
d3 −
rj
(∆t)3
dt3
(4.8)
δj (t) =
→
d4 −
rj
(∆t)4
4
dt
(4.9)
²j (t) =
→
d5 −
rj
(∆t)5
dt5
(4.10)
αj (t) =
e
74
→
Deste modo, o valor predito da posição para a j-ésima partı́cula, −
rj P (t + ∆t),
pode ser expresso por:
−
→
→
rj P (t + ∆t) = −
rj (t) + αj (t) + βj (t) + γj (t) + δj (t) + ²j (t)
(4.11)
Expandindo-se então os termos preditos, acima definidos, tem-se:
αjP (t + ∆t) = αj (t) + 2βj (t) + 3γj (t) + 4δj (t) + 5²j (t)
(4.12)
βjP (t + ∆t) = βj (t) + 3γj (t) + 6δj (t) + 10²j (t)
(4.13)
γjP (t + ∆t) = γj (t) + 4δj (t) + 20²j (t)
(4.14)
δjP (t + ∆t) = δj (t) + 5²j (t)
(4.15)
Uma vez que foram calculadas as posições preditas, procede-se à determinação
da força corrigida, de modo que esta pode ser reescrita da seguinte forma:
βjC (t + ∆t) =
1 Fj (t + ∆t)
(∆t)2
2!
mj
(4.16)
Logo, a diferença entre as equações 4.16 e 4.13 fornecerá a correção ∆βj (t+∆t)
do termo que contém a aceleração, de modo que:
βjC (t + ∆t) = βjP (t + ∆t) + ∆βj (t + ∆t)
(4.17)
Da mesma forma pode-se determinar os outros termos corrigidos, começando
pelas posições:
rjC (t + ∆t) = rjP (t + ∆t) +
3
∆βj (t + ∆t)
16
(4.18)
Determinando-se a seguir os termos αjC (t + ∆t), γjC (t + ∆t), δjC (t + ∆t) e
75
²C
j (t + ∆t), então:
αjC (t + ∆t) = αjP (t + ∆t) +
251
∆βj (t + ∆t)
360
(4.19)
γjC (t + ∆t) = γjP (t + ∆t) +
11
∆βj (t + ∆t)
18
(4.20)
1
δjC (t + ∆t) = δjP (t + ∆t) + ∆βj (t + ∆t)
6
(4.21)
1
∆βj (t + ∆t)
60
(4.22)
e
P
²C
j (t + ∆t) = ²j (t + ∆t) +
Os coeficientes
3 251 11 1
, , ,
16 360 18 6
e
1
60
são coeficientes ajustáveis visando, simulta-
neamente, acurácia e estabilidade para a solução. Os valores utilizados são tais que
representam o melhor compromisso entre a acurácia e estabilidade do algoritmo para
este tipo de problema. Algoritmos que se baseiem em outras ordens de aproximação
utilizam diferentes valores para os referidos coeficientes.
Resumidamente, constata-se que a idéia central de um algoritmo do tipo preditor − corretor é a comparação do valor da força calculada em duas estimativas
com diferentes ordens de precisão. Poder-se-ia fazer iterações sucessivas no algoritmo
preditor-corretor de forma a reduzir as correções dadas pela equação 4.17 e levar à
auto-consistência. No entanto, não há muita vantagem em tal procedimento uma vez
que o presente problema trata de sistemas fortemente correlacionados, ou seja, a cada
iteração uma nova força é calculada. Neste caso, resolveu-se empregar um intervalo
de tempo menor, pois mesmo que venha a se levar o algoritmo à consistência, não se
chega à trajetória exata do sistema, uma vez que o erro é de ordem n, num algoritmo
de ordem n. Pode-se também perceber que este algoritmo, pela utilização da terceira
lei de Newton, tem reduzido o tempo computacional à metade, uma vez que a força
que atua no átomo j, devido ao átomo i, é igual e oposta àquela atuante no átomo i,
devido ao átomo j.
76
4.2
4.2.1
Discussão e Resultados
Estudo do Campo Elétrico Médio
Nesta seção, será discutido preliminarmente o comportamento médio da in-
tensidade do campo elétrico gerado por superfı́cies fractais, restrigindo-se tal estudo
a sistemas com 100x100x100 pontos, cujos condutores base são superfı́cies geradas
mediante o algoritmo de DB ou por uma metodologia FBM. Esta discussão possui
duas motivações: a primeira é a necessidade de construção de um algoritmo que
determine o campo eletrostático em qualquer ponto do espaço e não apenas nos pontos
do “grid” utilizado para a solução da equação de Laplace. Deste modo, o estudo do
movimento de partı́culas carregadas não fica restrito a pontos do “grid”, ou seja, o
problema é extendido a uma situação contı́nua. O segundo aspecto está associado ao
estudo de como a rugosidade de um condutor base afeta o comportamento médio do
campo eletrostático gerado pelo mesmo.
Assim, para a determinação do módulo do campo eletrostático em qualquer
ponto do espaço, sabe-se que as componentes do vetor campo elétrico em cada ponto
do “grid”, nas direções dos eixos cartesianos representadas pelos vetores E~x , E~y e E~z ,
possuem módulos dados respectivamente por:
Ex =
{φ(i + 1, j, k) − φ(i, j, k)}
∆x
(4.23)
Ey =
{φ(i, j + 1, k) − φ(i, j, k)}
∆y
(4.24)
Ez =
{φ(i, j, k + 1) − φ(i, j, k)}
∆z
(4.25)
Nas equações 4.23, 4.24 e 4.25, tem-se que ∆x = ∆y = ∆z = 1, uma vez
que o espaçamento entre pontos vizinhos do “grid” foi considerado unitário. Neste
algoritmo, foi tida a preocupação de que o procedimento seja válido também para o
cálculo da intensidade do campo elétrico, nos pontos do “grid”, Egrid . Assim, estes
valores devem estar de acordo com a seguinte expressão:
77
1
Egrid = [(Ex )2 + (Ey )2 + (Ez )2 ] 2
(4.26)
Caso um determinado ponto P , de coordenadas (rx , ry , rz ), não pertença ao
“grid”, ou seja, esteja contido no interior de um cubo cujos vértices são pontos do
“grid”, a intensidade do campo elétrico é determinada considerando-se a influência
dos campos já calculados nestes oito pontos do vértice. Sejam as coordenadas que
localizam estes últimos, em relação à origem do sistema, denotadas por xn , yn e zn ,
com n = 0, 1 (ver Figura 4.1).
Figura 4.1: Esquema considerando o ponto P , cujo campo elétrico se quer determinar,
envolvido por oito pontos pertencentes ao “grid” regular.
Define-se também as quantidades rx0 , ry0 , rz0 , como:
rx0 = x0 − rx
(4.27)
ry0 = y0 − ry
(4.28)
rz0 = z0 − rz
(4.29)
Como o espaçamento, entre dois pontos vizinhos do “grid”, possui comprimento
78
considerado unitário, pode-se definir as quantidades rx1 , ry1 e rz1 de modo que:
rx1 = 1 − rx0
(4.30)
ry1 = 1 − ry0
(4.31)
rz1 = 1 − rz0
(4.32)
Portanto, denotando-se a distância entre o ponto cujo campo se quer determinar
e um dado ponto do vértice como rxn yn zn , com n = 0, 1, tem-se que:
q
rxn yn zn =
(rxn )2 + (ryn )2 + (rzn )2
(4.33)
A partir da discussão acima e das quantidades definidas até agora nesta seção,
a expressão para os módulos das componentes do campo elétrico, calculadas pelo
presente algoritmo, e denotadas por Exa = E1 , Eya = E2 e Eza = E3 são dadas por:
Eµ =
1
rnorm
"
1
X
Eµ(i+n,j+n,k+n)
n=0
rxn yn zn
#
Eµ(i+1,j,k+n) Eµ(i+n,j+1,k) Eµ(i,j+n,k+1)
+
+
+
rx1 y0 zn
rxn y1 z0
rx0 yn z1
(4.34)
Na expressão 4.34, onde µ = 1, 2, 3, os ı́ndices i, j e k localizam os pontos do
“grid” onde o campo elétrico foi determinado pelas expressões 4.23, 4.24 e 4.25 e o
fator de normalização, rnorm , é dado pela expressão:
rnorm =
1
X
n=0
"
1
rxn yn zn
+
1
rx1 y0 zn
+
1
rxn y1 z0
+
1
rx0 yn z1
#
(4.35)
Observa-se que a equação 4.34 leva em conta a dependência do campo em
um ponto qualquer com o inverso da distância deste ponto aos pontos do “grid” que
pertencem ao cubo e cujo campo elétrico já foi calculado a partir da solução numérica
da equação de Laplace. Um outro detalhe interessante é que esta equação determina
o valor do campo elétrico num ponto qualquer do espaço em função do campo nos
79
pontos do “grid”, para pontos no interior de um cubo definido pelos pontos do “grid”.
Por fim, o presente algoritmo, atende à preocupação já mencionada para o cálculo do
campo no caso particular em que as coordenadas correspondem às de um dos pontos
do “grid”, ou seja:
1
1
[(Exa )2 + (Eya )2 + (Eza )2 ] 2 = [(Ex )2 + (Ey )2 + (Ez )2 ] 2
(4.36)
Uma vez que o campo elétrico é determinado em um ponto qualquer do espaço,
resolveu-se investigar como o módulo do campo elétrico médio, associado a cada
superfı́cie equipotencial φi , se comporta com respeito à distância média desta última
ao condutor base. Esta distância média foi definida no Capı́tulo 3, a partir da equação
3.30. Este estudo foi motivado por se buscar interpretar como condutores base, com
diferentes rugosidades, afetam o comportamento médio do campo eletrostático. Foram
consideradas as superfı́cies condutoras geradas a partir dos algoritmos DB e FBM com
tamanhos 100x100. A intensidade do campo elétrico médio, < E >(φi ) , associada a
cada superfı́cie equipotencial é dada por:
L2
1
1 X
< E >(φi ) = 2
[(Emx )2 + (Emy )2 + (Emz )2 ] 2
L m=1
(4.37)
Na equação 4.37 o ı́ndice m indica um ponto pertencente a uma dada superfı́cie
equipotencial, cujo campo elétrico médio se quer determinar. Como o presente estudo
foi restrito a sistemas de dimensões 100x100, para cada superfı́cie equipotencial,
serão considerados, para o cálculo da média, 10000 pontos. A figura 4.2 mostra o
comportamento do desvio quadrático médio como função da escala, para as superfı́cies
condutoras base utilizadas como fronteiras para solução numérica da equação de
Laplace.
Observa-se que se tratam de fractais, de modo que a superfı́cie condutora
gerada por DB, BAL100(1), apresenta dimensão fractal Df 1 = (2, 482 ± 0, 005). As
superfı́cies condutoras base geradas pela metodologia FBM, denotadas por F BM 100(2)
e F BM 100(3), apresentam dimensões fractais Df 2 = (2, 486 ± 0, 007) e Df 3 =
(2, 408 ± 0, 007) respectivamente (ver Figura 4.3). Para este estudo o método RMS
80
Figura 4.2: Comportamento do desvio quadrático médio como função da escala, em
escala logarı́tmica, tı́pico de fractais auto-afins. Os condutores base foram gerados por
DB e pela metodologia FBM. Observa-se que a rugosidade é maior para o condutor
gerado por deposição balı́stica para todas as escalas consideradas. Para as superfı́cies
geradas pela metodologia FBM as rugosidades são diferentes entre si, de modo que
considerou-se uma superfı́cie de rugosidade intermediária e uma de menor rugosidade
comparativamente.
foi considerado para determinação dos expoentes de rugosidade dada à sua maior
sensibilidade, conforme já mencionado anteriormente. A rugosidade do condutor base
gerado por DB, BAL100(1), foi 3, 38. Para os condutores gerados pela metodologia
FBM, F BM 100(2) e F BM 100(3), as rugosidades foram 3, 38 e 1, 44.
A figura 4.4 mostra o comportamento da intensidade do campo elétrico médio
como função da distância média, considerando-se os condutores base já descritos.
Trata-se de um resultado a priori não esperado, uma vez que verifica-se um crescimento inicial do campo elétrico médio com o aumento da distância média. A influência
da rugosidade neste comportamento pode ser percebida no intervalo dos valores associados ao campo médio até que a sua magnitude se torne constante. Neste caso,
este intervalo é claramente maior considerando-se o condutor base gerado por DB e
menor considerando-se o gerado por uma metodologia FBM com a menor rugosidade,
ou seja, 1, 02. Outro dado a observar neste resultado é a existência de um campo
médio máximo correspondente à uma dada distância média.
No gráfico 4.5, o comportamento da distância média das superfı́cies equipo-
81
Figura 4.3:
Superfı́cies base condutoras utilizadas para investigação do
comportamento do campo eletrostático médio. A região representada em branco
representa protuberâncias com maiores alturas, enquanto que a região enegrecida
representa pontos de menor altura.
Figura 4.4: Comportamento da intensidade do campo elétrico médio como função
da distância média. O cálculo do campo médio bem como da distância média ao
condutor base, foi feito considerando-se as superfı́cies equipotenciais cujos potenciais
elétricos são 99, 98, 97, 95, 90, 90 - n (n=10,20,..,60), 25, 15. Está indicado também
neste gráfico a rugosidade, W (L), de cada condutor base considerado neste estudo e
já mencionada anteriormente.
82
tenciais ao condutor base é mostrado para os condutores BAL100(1), FBM100(2) e
FBM100(3), indicando para que valores de distâncias médias, o campo médio assume
valor máximo. Coincidentemente, o potencial associado a este valor máximo é o
mesmo para os três condutores base, ou seja, φ = 70.
Figura 4.5: Comportamento da distância média de cada superfı́cie equipotencial ao
correpondente condutor base, como função do potencial elétrico. Para este cálculo
foram consideradas as superfı́cies equipotenciais cujos potenciais elétricos são 99, 98,
97, 95, 90, 90 - n (n=10,20,..,60), 25, 15. Está indicado também neste gráfico os
potenciais elétricos, cujo campo elétrico médio assume a máxima intensidade, bem
como os correspondentes valores. Observa-se uma coincidência na equipotencial,
cujo campo elétrico médio apresenta máximo valor, no caso, φi = 70, para os três
condutores base considerados.
Os resultados indicam que a intensidade do campo médio máximo possui
maior magnitude considerando-se superfı́cies base mais rugosas, sugerindo-se portanto
que regiões onde o campo elétrico possui maior intensidade no condutor base, como
conseqüência do efeito das pontas, afetam a topologia das superfı́cies equipotenciais
distantes trazendo influências relevantes no cálculo do campo médio. É observado,
mais claramente para o caso do condutor base gerado por DB, que após o valor máximo
assumido pelo campo elétrico médio, este decai suavemente até que se torne uniforme.
Para os outros condutores base, este efeito é também verificado, apresentando contudo
um decaimento ainda mais suave, comparado com o caso da superfı́cie base gerada
83
por DB.
Outra questão que é possı́vel de ser observada pelo exame da Figura 4.4, é o
fato de que a intensidade do campo médio para uma superfı́cie equipotencial próxima
ao condutor base, por exemplo φi = 99, apresenta-se maior para o condutor menos
rugoso. Este resultado sugere que, quanto mais rugoso seja um condutor base, maior
a possibilidade, no cálculo do campo elétrico num ponto qualquer do espaço (equação
4.34), da soma para cada componente (x,y e z) possuir termos que levem a soma a
apresentar menores valores. Portanto, os resultados indicam que, considerando-se um
condutor base mais rugoso, associa-se à uma superfı́cie equipotencial próxima, um
campo elétrico médio cuja intensidade é relativamente menor.
4.2.2
Movimento de Partı́culas não Interagentes em um Campo Elétrico Irregular
Neste trabalho, é interessante discutir também a dinâmica de partı́culas carre-
gadas sujeitas a um campo elétrico irregular, produzido por condutores que apresentem
geometria fractal. Estas partı́culas foram consideradas como não interagentes entre si,
sendo esta opção motivada pelo fato de que em emissores mais comuns, tipo “spindt”,
por exemplo [39], a densidade das partı́culas emitidas (elétrons ou ı́ons) é geralmente
baixa o suficiente para não serem incluı́dos efeitos de carga espacial. Buscou-se
então analisar que possı́veis informações, a respeito do campo eletrostático local bem
como da geometria do condutor base, podem ser obtidas a partir de uma análise
das trajetórias descritas e das grandezas fı́sicas dinâmicas associadas às partı́culas,
como por exemplo a energia cinética. Para isto foi adotado o método PreditorCorretor, já descrito anteriormente, aplicado a uma dinâmica clássica de partı́culas.
A utilização deste método é justificada, uma vez que o sistema em estudo é fortemente
correlacionado, ou seja, para cada passo de tempo que se considere, cada partı́cula
estará sujeita a uma força eletrostática diferente. Considerou-se neste estudo que a
velocidade inicial associada a cada partı́cula possui valor nulo, uma vez que, assim
sendo, as trajetórias descritas apresentam configurações que se assemelham às linhas
84
de força que caracterizam o campo eletrostático estudado (ver Figura 4.6).
Figura 4.6: (a) Partı́culas abandonadas na imediata vizinhança de uma protuberância
paraboloidal. (b) Trajetórias das partı́culas calculadas pelo método de dinâmica
molecular (preditor-corretor).
Uma conseqüência direta deste fato é que a forma da trajetória descrita por
uma dada partı́cula praticamente independe da razão entre a carga e a sua massa
correspondente, resultado este que pode ser verificado a partir da Figura 4.7.
Para avaliar como as irregularidades locais associadas ao condutor base podem
afetar o comportamento relativo à dinâmica de partı́culas, a energia cinética de cada
partı́cula foi calculada como função da posição correspondente à coordenada z, para
cada passo de tempo considerado. Considerou-se neste cálculo, para o condutor base
gerado por DB, que as coordenadas que definem uma superfı́cie equipotencial próxima
ao mesmo, no caso φ = 99, fossem correspondentes às posições iniciais associadas às
partı́culas. Estas últimas totalizavam 10000, recobrindo portanto a referida superfı́cie
base. Considerou-se para este cálculo que o plano de chegada destas partı́culas era
tal que z = 40. A partir da Figura 4.8, observa-se tal comportamento para algumas
partı́culas, considerando-se como condutor base, BAL100(1). Nota-se que àquelas
partı́culas que são abandonadas de regiões que possuem valores de z acima de 32, ou
seja, mais protuberantes, estão associadas uma diminuição no gradiente da energia
cinética com respeito à direção que define o eixo z, com a evolução temporal. Isto é
justificado pelo fato de se considerar o efeito das pontas, ou seja, tais partı́culas, a
85
Figura 4.7: À esquerda são mostradas as trajetórias, de uma dada partı́cula
considerando-se as razões entre as cargas e as massas distintas, com a superfı́cie
base Bal100(1) mostrada na Figura 4.3. Verifica-se que para diferentes razões mq ,
as trajetórias praticamente se sobrepõem. À direita, essas trajetórias são mostradas
para o caso em que a superfı́cie base é gerada por uma metodologia FBM, FBM100(2).
Nos dois casos, nos planos x=0, y=0 e z=0 são mostradas as projeções correspondentes
da trajetória descrita.
princı́pio, estão submetidas a campos mais intensos que os associados às regiões de
depressão. Para partı́culas abandonadas de posições iniciais associadas a menores
valores de z observa-se uma inversão no comportamento do gradiente da energia
cinética com relação a z, justificando-se portanto a presença de campos relativamente
menores, associados a regiões menos protuberantes. Este efeito de inversão no entanto,
também foi verificado considerando-se como condutores base, os gerados por uma
metodologia FBM (ver Figura 4.9), sendo contudo menos acentuado. Considerou-se
neste caso o plano de chegada z = 30. Este resultado sugere que tal efeito é verificável
mais freqüentemente para superfı́cies base menos suaves. Vale ressaltar ainda que o
efeito de distância ao plano de chegada, faz com que às partı́culas cujas posições
iniciais estejam associadas às regiões mais protuberantes, estejam associadas menores
magnitudes de energia cinética, na chegada. Esta hipótese de descorrelação é evitada
naturalmente ao se considerar as posições iniciais associadas às 10000 partı́culas como
sendo um plano horizontal inicial, zi , localizado imediatamente acima do condutor
base, ou seja, zi =f (xgrid , ygrid ), onde xgrid e ygrid são coordenadas pertencentes ao
“grid” regular. Por fim, calcula-se a energia cinética, Ecf , das partı́culas em um
86
plano de chegada, zf >zi , não muito afastado do plano inicial. Assim, a “superfı́cie
de energia” Ecf = g(xgrid , ygrid ) é construı́da, refletindo a conservação da energia
mecânica. Naturalmente, neste caso, o efeito de descorrelação entre a energia cinética
e o campo elétrico é eliminado uma vez que a componente do vetor posição na direção
k̂ é o mesmo para todas as partı́culas que chegam ao plano de chegada. Na Figura
4.10, mostra-se uma comparação visual entre a superfı́cie de energia definida acima,
e a distribuição dos potenciais elétricos no plano que define as posições iniciais das
partı́culas zi = f (xgrid , ygrid ), utilizando-se como superfı́cie base BAL100(1).
Figura 4.8: Comportamento da energia cinética de algumas partı́culas representadas
por P − 2, P − 2319, ..., como função de z considerando-se como superfı́cie base,
BAL100(1). Com relação ao gráfico amarelo, verificou-se que a partı́cula associada
antes de chegar ao plano z = 40 ultrapassa lateralmente a região de integração
numérica.
Para verificar como a energia cinética das partı́culas, no “plano de chegada”
pode trazer informações a respeito da geometria do condutor base, utilizou-se o
método RMS. Para este estudo, considerou-se que um número de partı́culas cujas
equações do movimento serão resolvidas é equivalente ao tamanho do sistema, ou
seja, 10000. Estas partı́culas, foram abandonadas do repouso possuindo coordenadas
iniciais iguais às coordenadas que definem a superfı́cie equipotencial φ = 99. Calculouse a energia cinética, Ecf , das partı́culas em um plano de chegada, zf , não muito
87
Figura 4.9: Comportamento da energia cinética de algumas partı́culas representadas
por P − 2, P − 2319, ..., como função de z considerando-se como processo de geração
do condutor base a metodologia FBM. Neste caso o plano de chegada foi z = 30 e a
superfı́cie base, FBM100(2).
Figura 4.10: Na Figura à esquerda é mostrada a superfı́cie de energia Ecf =
g(xgrid , ygrid ). Os pontos mais claros representam maiores magnitudes de energia
cinética, enquanto que os mais escuros, o oposto. À direita a distribuição dos
potenciais no plano zi = f (xgrid , ygrid ) é mostrado. Foi utilizado neste cálculo o
condutor base BAL100(1) mostrado na Figura 4.3.
88
afastado do condutor base. Assim, uma “superfı́cie de energia” Ecf = g(xgrid , ygrid ),
onde xgrid e ygrid são coordenadas pertencentes ao “grid” regular, foi construı́da
objetivando-se determinar como se dá a distribuição da energia cinética no plano
de chegada. Os “planos de chegada” como já exposto, foram considerados zf = 40
para o condutor base BAL100(1) e zf = 30 para os condutores base F BM 100(2)
e F BM 100(3). Naturalmente, as energias cinéticas associadas a cada partı́cula não
serão iguais no plano de chegada considerado, uma vez que este último não é uma
superfı́cie equipotencial. Conseqüentemente, este plano possui uma distribuição de
potenciais elétricos que guarda informações a respeito da geometria do condutor base.
Na figura 4.11 tem-se as “superfı́cies de energia” considerando-se os condutores base
já descritos. Eliminou-se a anticorrelação existente devido aos efeitos de distância
entre o campo elétrico inicial e a energia cinética no plano de chegada, a fim de se
facilitar a identificação de regiões que compõem a imagem da “superfı́cie de energia”
que apresentam semelhanças geométricas com o condutor base. Na figura 4.12 são
mostradas as superfı́cies equipotenciais calculadas, φ = 99, considerando-se os condutores base mostrados na Figura 4.3.
Na Figura 4.13 observa-se o comportamento do desvio quadrático médio (RMS)
normalizado como função da escala para as superfı́cies equipotenciais φ = 99 e
para as correspondentes “superfı́cies de energia” apresentadas nas Figuras 4.11 e
4.12 respectivamente. Como se pode verificar, este resultado sugere a existência
de comprimentos de de escala para os quais os valores de desvio quadrático médio
normalizados são praticamente coincidentes para as superfı́cies equipotenciais φ = 99
e para as correspondentes “superfı́cies de energia”. Este aspecto justifica a semelhança
apresentada entre as figuras 4.11 e 4.12 correspondentes, confirmando a hipótese
de que a medida da energia cinética de partı́culas em um campo irregular pode
trazer informações geométricas que, a depender do comprimento de escala utilizado, podem ser relevantes para se determinar propriedades de rugosidade do condutor
base estudado. Outro aspecto interessante a se verificar é que as inclinações, α, no
ajuste linear calculado, considerando-se como escala limite o valor 20, são menores
para as “superfı́cies de energia” implicando-se em uma maior dimensão fractal. No
89
Figura 4.11: “Superfı́cies de energia”, Ecf = g(xgrid , ygrid ), considerando-se o
movimento de partı́culas carregadas não interagentes num campo elétrico gerado
pelos condutores base BAL100(1), FBM100(2) e FBM100(3).
As partı́culas
foram abandonadas do repouso com posições iniciais definidas pelas coordenadas,
(xgrid , ygrid , zφ=99 ), das correspondentes superfı́cies equipotenciais φ = 99. Os pontos
mais claros representam maiores magnitudes de energia cinética, enquanto que os mais
escuros, o oposto.
90
Figura 4.12: Superfı́cies equipotenciais, φ = 99, calculadas a partir da solução
numérica da equação de Laplace, para os condutores base BAL100(1), FBM100(2) e
FBM100(3). Os pontos mais claros representam regiões mais protuberantes, enquanto
que os mais escuros, o oposto.
91
entanto, considerando-se para a determinação da dimensão fractal valores iniciais e
finais de escala maiores, as inclinações apresentam valores mais próximos, sugerindo
dimensões fractais mais próximas. Este aspecto é interessante, pois permite a obtenção
da dimensão fractal do condutor base, a partir da distribuição das energias cinéticas
das partı́culas, considerando-se como pontos de ajuste no método RMS os que possuem
escalas maiores e que não pertençam aos pontos que compõem a saturação da curva.
Uma análise adicional foi feita para uma comparação quantitativa das superfı́cies equipotenciais bem como das correspondentes “superfı́cies de energia” utilizandose o método da análise de Fourier discutida no Capı́tulo 2. Nas Figuras 4.14, 4.15 e
4.16 são mostrados os comportamentos do quadrado das amplitudes como função das
freqüências, constituindo-se o espectro de potência, para as superfı́cies ENBAL100(1)
e BAL100(1) - φ = 99, ENFBM100(2) e FBM100(2) - φ = 99 e ENFBM100(3)
e FBM100(3) - φ = 99 respectivamente. São mostrados também o ajuste linear
utilizando-se o método dos mı́nimos quadrados e os valores obtidos para os coeficientes
lineares, A, e angulares, β.
É interessante observar que o processo de ajuste linear de muitos pontos no
espectro de Fourier, através do método dos mı́nimos quadrados, inclui poucos pontos
que se afastam significativamente da tendência principal do ajuste.
Contudo, a
presença de picos maiores na transformada de Fourier é visualmente evidente nas
Figuras 4.14, 4.15 e 4.16 e ignorada no ajuste linear. Isto significa que cada pico
corresponde a uma única senóide sobreposta dominando a elevação da superfı́cie. No
entanto, esta é a menor componente do gráfico log(Amplitude)2 vs. log(F reqüência).
Outro aspecto é que a inclinação da reta de ajuste, é fortemente afetada pelo deslocamento dos pontos no fim do conjunto de dados, ou seja, as flutuações correspondentes
às altas freqüências podem levar a um crescimento nos valores da direita do gráfico
log(Amplitude)2 vs. log(F reqüência), levando a um deslocamento dos pontos correspondentes para cima. Este fato reduz a inclinação da reta de ajuste.
Um problema mais sutil e talvez o mais comum durante o ajuste, tem a ver com
os pontos da esquerda no espectro de potência [25]. Os termos de baixa freqüência
estão em menor número e conseqüentemente com maior espaçamento em relação ao
92
Figura 4.13: (a) Comportamento do desvio quadrático médio como função da escala,
considerando-se a superfı́cie de energia ENBAL100(1) e a superfı́cie equipotencial
BAL100(1) - φ = 99 mostradas nas figuras 4.11 e 4.12. (b) Comportamento do
desvio quadrático médio como função da escala, considerando-se a superfı́cie de energia
ENFBM100(2) e a superfı́cie equipotencial FBM100(2) - φ = 99 mostradas nas figuras
4.11 e 4.12. (c) Comportamento do desvio quadrático médio como função da escala,
considerando-se a superfı́cie de energia ENFBM100(3) e a superfı́cie equipotencial
FBM100(3) - φ = 99 mostradas nas figuras 4.11 e 4.12. A circunferência que contém
alguns dados indica a presença de interseções.
93
Figura 4.14: Espectros de potência para as superfı́cies BAL100(1) - φ = 99 e
ENBAL100(1), representados pelas cores preta e vermelha respectivamente. Os
coeficientes angular e linear são também explicitados no ajuste linear, utilizando-se o
método dos mı́nimos quadrados.
Figura 4.15: Espectros de potência para as superfı́cies FBM100(2) - φ = 99 e
ENFBM100(2), representados pelas cores preta e vermelha respectivamente. Os
coeficientes angular e linear são também explicitados no ajuste linear, utilizando-se o
método dos mı́nimos quadrados.
94
Figura 4.16: Espectros de potência para as superfı́cies FBM100(3) - φ = 99 e
ENFBM100(3), representados pelas cores preta e vermelha respectivamente. Os
coeficientes angular e linear são também explicitados no ajuste linear, utilizando-se o
método dos mı́nimos quadrados.
eixo das freqüências em escala logarı́tmica, controlando ainda mais a inclinação da
reta de ajuste. Contudo, estes termos incluem provavelmente a informação da forma
total da superfı́cie associada ao processo de construção. O primeiro termo da série de
Fourier é simplesmente a altura média da superfı́cie, que sendo arbitrária depende da
escolha de um ponto de referência.
Apesar destas limitações, pode-se verificar que, no ajuste linear, os módulos
das inclinações obtidas com relação às superfı́cies de energia sofrem uma redução em
relação às correspondentes superfı́cies equipotenciais φ = 99. Este fato sugere que as
superfı́cies de energia são mais irregulares uma vez que a inclinação no ajuste linear
está associada ao cálculo da dimensão fractal, Df pela expressão já apresentada no
Capı́tulo 2:
Df =
(6 + β)
2
(4.38)
Observa-se que as inclinações, β, associadas às superfı́cies equipotenciais BAL95
100(1) − φ = 99, F BM 100(2) − φ = 99 e F BM 100(3) − φ = 99, são tais que
βBAL100(1)−φ=99 > βF BM 100(2)−φ=99 > βF BM 100(3)−φ=99 . Esta desigualdade é também
evidente ao se comparar as inclinações obtidas para as correspondentes superfı́cies de
energia. Este fato confirma a correlação existente entre as geometrias das superfı́cies
de energia e as geometrias dos condutores base correspondentes.
Por fim, neste trabalho, é discutido também como se dá a distribuição das
partı́culas em um plano de chegada distante do condutor base, que neste caso corresponde ao plano zf = 100.
Este plano define o limite superior na região de
integração onde o potencial elétrico foi definido como condição de fronteira, φ = 0. As
posições iniciais associadas às 10000 partı́culas que foram abandonadas do repouso,
correspondem às coordenadas que definem a superfı́cie equipotencial φ = 99. Foram
determinadas as densidades superficiais de partı́culas, σ(Ei ), ou seja, o número de
partı́culas por unidade de área do “grid” no plano de chegada. O seu comportamento
foi analisado como função do campo elétrico associado à posição inicial de cada
partı́cula, Ei . Na Figuras 4.17 e 4.18, são mostradas a superfı́cie σ(Ei )=h(xgrid , ygrid )
e o comportamento de σ(Ei ) como função do campo elétrico Ei para as superfı́cies
BAL100(1) e FBM100(2).
Fazendo-se uma comparação visual entre as Figuras 4.17 e 4.18 (a) e as correspondentes superfı́cies equipotenciais indicadas na Figura 4.12 verifica-se uma anticorrelação. Às regiões associadas ao condutor base que possuem um campo elétrico
mais intenso, associa-se nas figuras 4.17 e 4.18 (a) regiões cuja densidade superficial de
partı́culas é menor. Estas últimas são representadas pela coloração negra. Este fato
está associado à presença de maiores intensidades nas componentes Ex e Ey do campo
elétrico em regiões mais protuberantes do condutor base. Isto faz com que a trajetória
correspondente a cada partı́cula seja mais intensamente desviada da direção vertical.
Observa-se claramente que existem regiões onde a densidade superficial de partı́culas
é maior, estando relacionada a campos elétricos menos intensos gerados pelo condutor
base. Este aspecto sugere que o grau de irregularidade apresentada por um condutor
base pode ser importante para que, a partir da dinâmica de elétrons, procure-se uma
possı́vel otimização de dispositivos emissores em escalas nanométricas. Esta vantagem
96
Figura 4.17: (a) Superfı́cie σ(Ei )=h(xgrid , ygrid ) considerando-se o condutor base
BAL100(1). As regiões brancas representam uma maior densidade superficial de
partı́culas no plano z = 100. (b) Comportamento da densidade superficial de
partı́culas como função da intensidade do campo elétrico associado à posição inicial
de cada partı́cula. .
Figura 4.18: (a) Superfı́cie σ(Ei )=h(xgrid , ygrid ) considerando-se o condutor base
FBM100(2). As regiões brancas representam uma maior densidade superficial de
partı́culas no plano z = 100. (b) Comportamento da densidade superficial de
partı́culas como função da intensidade do campo elétrico associado à posição inicial
de cada partı́cula.
97
se refere à captura de poucos elétrons, tendo-se portanto correntes elétricas de baixa
intensidade permitindo assim diminuições consideráveis nas dissipações de energia.
Um fato curioso que pode ser observado nas Figuras 4.17 e 4.18 (b) está associado às
intensidades dos campos elétricos nos pontos que definem a superfı́cie equipotencial
φ = 99 para os condutores base BAL100(1) e F BM 100(2). Apesar da superfı́cie
F BM 100(2) ser mais homogênea e menos rugosa em relação à superfı́cie BAL100(1),
verifica-se a presença de campos mais intensos para o condutor F BM 100(2), resultado
este de acordo com os apresentados no estudo do campo elétrico médio.
98
Capı́tulo 5
Conclusões e Perspectivas
Neste trabalho estudou-se o campo eletrostático gerado por condutores que
apresentam irregularidades com uma geometria fractal. No estudo para sistemas com
1D+1 dimensões, a equação de Laplace foi resolvida numericamente numa região entre
um perfil condutor base com geometria auto-afim e construção iterativa bem definida,
e uma linha distante. As condições de Dirichlet foram atribuı́das mantendo-se entre
o condutor base rugoso e o distante, retilı́neo, uma diferença de potencial constante.
Posto que não é possı́vel uma solução analı́tica deste problema, os perfis equipotenciais
foram determinados por interpolação linear e as irregularidades associadas foram
representadas pelos expoentes de rugosidade, estes últimos determinados pelo método
do Semivariograma. Procurou-se, nesta parte do trabalho, avaliar como a dimensão
fractal se comporta em função do potencial elétrico para sistemas de tamanhos fixos.
Investigou-se também como a mudança do tamanho de um sistema afeta a dimensão
fractal dos perfis equipotenciais. Os resultados indicam que, para um tamanho fixo, a
dimensão fractal varia com o potencial elétrico seguindo uma relação que não obedece
leis exponenciais ou leis de potência. Um aumento da distância média da linha
equipotencial ao condutor base acarreta um aumento dos expoentes de rugosidade
e, conseqüentemente, uma diminuição das correspondentes dimensões fractais. Com
relação aos efeitos de influência do tamanho do sistema na dimensão fractal de uma
dada linha equipotencial φi , os resultados indicam que o aumento do tamanho do
sistema acarreta um aumento da dimensão fractal correspondente. Este fato sugere
99
que, no limite de sistemas infinitos, a dimensão fractal depende da distância média
de cada linha equipotencial ao condutor base. Para sistemas infinitos, todas as linhas
equipotenciais apresentariam a mesma dimensão fractal do condutor base considerado.
Outra questão de interesse, ainda para sistemas com 1D+1 dimensões, é o
comportamento da distância média de cada linha equipotencial ao condutor base
como função do potencial elétrico. Foram verificados comportamentos não lineares
mais acentuados para potenciais intermediários, confirmando-se os resultados obtidos
em trabalhos anteriores (Dias Filho et al [32]), onde perfis randômicos base foram
utilizados.
Neste trabalho foi também abordada a solução numérica da equação de Laplace
para a determinação das superfı́cies equipotenciais, considerando-se sistemas com
2D+1 dimensões. O método de Liebmann foi portanto extendido para 2D+1 dimensões, ou seja, a região de integração é o espaço tri-dimensional. Foram considerados como condutores base, neste estudo, uma superfı́cie resultado de uma deposição
balı́stica e uma superfı́cie gerada por uma metodologia FBM. No processo de geração
das superfı́cies por deposição balı́stica foi necessário, para a obtenção da superfı́cie
resultante, se considerar algumas aproximações particulares. Uma está associada ao
processo de construção devido à correlação lateral existente no processo de deposição.
A outra está associada à saturação da rugosidade.
Para quantificação das irregularidades das superfı́cies equipotenciais foram
utilizados os métodos RMS e Semivariograma. Os resultados permitem concluir que
o método RMS é mais sensı́vel que o método do Semivariograma. Os domı́nios
de validade destes métodos foram investigados, determinando-se um expoente de
rugosidade médio como limite, em torno do qual os métodos perdem sua validade
de aplicação. Comparando-se o comportamento da dimensão fractal de superfı́cies
equipotenciais com o respectivo potencial elétrico, para os condutores base discutidos
neste estudo, todos com o mesmo tamanho, os resultados indicam um decréscimo mais
acentuado da dimensão fractal com o decréscimo do potencial elétrico para condutores
base gerados por uma metodologia FBM. Este resultado sugere que, sendo o condutor
base FBM mais suave, as superfı́cies equipotenciais perdem mais rapidamente a
100
correlação geométrica com o condutor base resultando em um campo eletrostático
praticamente uniforme.
Um comportamento não trivial entre a distância média e o potencial elétrico
foi encontrado, no limite em que < d(φ) >→0. Esta transição da não linearidade,
para situações onde o comportamento de < d(φ) > com φ é mais aproximadamente
linear, sugere que as variações das rugosidades das superfı́cies equipotenciais próximas
ao condutor base com o potencial elétrico são mais acentuadas enquanto que para
potenciais elétricos intermediários estas variações são mais suaves.
Analisou-se também a dinâmica de partı́culas carregadas não interagentes,
sujeitas a campos elétricos com significativa variação local. Utilizou-se o método
preditor-corretor para a determinação das quantidades dinâmicas. Para a simulação
computacional houve a necessidade de construção de um algoritmo que possibilitasse
determinar o campo eletrostático em qualquer ponto do espaço, o que leva a um
tratamento do problema de modo equivalente ao de um contı́nuo. Isto possibilitou
que fosse investigado o comportamento médio do campo elétrico, associado a cada
superfı́cie equipotencial, como função da distância média de cada equipotencial ao
condutor base. Foram então considerados condutores base com rugosidades diferentes.
Os resultados indicam a existência de um campo elétrico médio máximo, que apresenta
dependência com a rugosidade do condutor base. Para condutores mais rugosos
campos elétricos médios máximos de maior intensidade relativa, foram assim encontrados. Os resultados indicam também que a medida de grandezas dinâmicas,
associadas às partı́culas carregadas não interagentes, como a energia cinética, pode
conter/revelar, a depender do comprimento de escala considerado, informações a
respeito das irregularidades do condutor base. Os métodos RMS e a análise de Fourier
foram utilizados para quantificar estas informações que, para a metodologia utilizada
neste trabalho, consistiu na construção de “superfı́cies de energia”.
Este conjunto de resultados, além de representar uma contribuição cientı́fica
para o entendimento de propriedades do campo elétrico gerado por condutores de
geometria bastante irregular, também tem o objetivo de contribuir na investigação
e interpretação do comportamento de dispositivos e de processos que ocorrem em
101
superfı́cies irregulares. Isto ocorre em diversas áreas da Fı́sica da Matéria Condensada,
como no estudo do processo de adsorção de partı́culas em superfı́cies irregulares,
no processo de emissão por campo associado a nanoestruturas irregularidades, os
processos de transporte de elétrons em nanoestruturas metálicas, etc... . Como
perspectiva futura, um estudo do transporte de nanoestruturas em campos elétricos
irregulares poderá ser feito, objetivando-se investigar durante a acomodação de tais
estruturas sobre o condutor base, efeitos de distribuição eletrônica, visando a otimização e a construção de novos materiais.
Outro estudo que merece atenção, como perspectiva de aplicação futura, é o
estudo da emissão de elétrons a partir de condutores irregulares. Considerando a
variação do campo elétrico local, é possı́vel, mediante a equação de Fowler-Nordheim
[40], avaliar a área de extração do emissor. O presente estudo permite determinar
fatores de amplificação do campo elétrico em função da rugosidade do condutor base
considerado. Esta análise permite verificar a hipótese de que propriedades associadas
a emissão por campo de um grande número de materiais estariam conectadas com
a estrutura fractal das superfı́cies. Em caso de confirmação desta hipótese, o fator
de amplificação deveria ser considerado para o cálculo da densidade de corrente de
emissão.
102
Referências Bibliográficas
[1] M. V. Eliseev, O. B. Isayeva, A. G. Rozhnev, N. M. Ryskin, Elsev. Sci. 45,
871-877 (2001).
[2] S. R. Salinas, Introdução à Fı́sica Estatı́stica (EDUSP, 1997).
[3] J. Feder Fractals (Plenum Press, New York, 1988).
[4] A. L. Barabási, H. E. Stanley, Fractal Concepts in Surface Growth (Cambridge
University Press, 1995).
[5] S. S. Brenner, J. T. McKinney, Surf. Sci. 23, 88 (1970).
[6] Wojciechowski K F Europhys. Lett. 38, 135 (1997).
[7] B. B. Mandelbrot, The Fractal Geometry of Nature (New York, W. H. Freeman
and Company, 1975).
[8] J. W. Dauben, Georg Cantor: His Mathematics and Philosophy of the Infinite.
(Cambridge, Mass, 1979)
[9] H. von. Koch, Sur une courbe continue sans tangente, obtenue par une
construction gomtrique lmentaire. (Archiv fr Matemat., Astron. och Fys. 1, 681702, 1904).
[10] F. Barnsley, J. S. Geronimo, and A. N. Harrington, Geometrical and electrical
properties of some Julia sets (Lecture Notes in Pure and Applied Mathematics,
Vol. 92, 1988).
[11] E. Lorenz, J. Atmos. Sci 20, 130-141 (1963).
103
[12] T. Y. Li, J. A. Yorke, Am. Math. Monthly 82, 985-992 (1975).
[13] R. M. May, J. Theor. Biol., 51, 511524 (1975).
[14] F. C. Hoppensteadt, J.M. Hyman, (Courant Institute, New York University:
preprint, 1975).
[15] B. B. Mandelbrot. Fractals - Form, Chance, Dimension. (San Francisco,
California, 1977).
[16] K. Falconer, Fractal Geometry - Mathematical Foundation and Applications
(England, 1990).
[17] R. R. Goldberg, Methods of Real Analysis (Blaisdell, New York,1963).
[18] E. L. Lima, Análise Real (V1. IMPA, Rio de Janeiro, 1989).
[19] W. Sierpinski, Elem. Math., 15, 73-74 (1960).
[20] B. B. Mandelbrot. Objetos Fractais (Lisboa, 1991).
[21] B. Dubuc, J. F. Quiniou, C. Roques, C. Tricot, Phys. Rev. A., 39, 15001512
(1989).
[22] R. F. S. Andrade, D. O. Cajueiro, C. S. Ferreira, Physica A, 295, 323-332 (2001).
[23] D. Avnir, The Fractal Approach to Heterogeneous Chemistry (Hebrew University
of Jerusalem, 1989).
[24] J. G. V. Miranda, Análisis Fractal Del Microrrelieve del Suelo. (La Coruña:
Espanha, 2000. 258 p. Tese (Doutorado)).
[25] J. C. Russ, Fractal Surfaces (New York and London Press, 1994).
[26] J. FEDER, Fractals (New York Plenum Press, 1988).
[27] G. S. P. Miller, The Definition and Rendering of Terrain Maps, ACM
SIGGRAPH Computer Graphics, 20, 39-48 (1986).
104
[28] P. Meakin, Fractals, Scaling and Growth Far from Equilibrium (Cambridge
University Press, New York, 1998).
[29] C. F. Gerald, P. O. Wheatley, Applied Numerical Analisys, Third Edition,
Addison-Wesley, (New York, 1985).
[30] T. T. Tsong, Atom-Probe Field Ion Microscopy (Cambridge University Press,
1990).
[31] D. O. Cajueiro, V. A. de A. Sampaio, C. M. C. de Castilho e R. F. S. Andrade,
J. Phys.: Condens. Matter 11 4985 (1999).
[32] H. de O. D. Filho, C. M. C. de Castilho, J. G. V.Miranda, R. F. S. Andrade,
Physica A , 342, 388 (2004).
[33] T. A. de Assis, F. de B. Mota, J. G. V. Miranda, R. F. S. Andrade, H. de O.
Dias Filho e C. M. C. de Castilho, J. Phys.: Condens. Matter 18 3393 (2006).
[34] P. Woodruff and T. A. Delchar 1994 Modern Techiniques os Surface Science 2nd
ed (Cambridge: Cambridge University Press).
[35] I. Brodie Phys. Rev. B,51, 13660 (1995).
[36] K. F. Wojciechowski Europhys. Lett,38, 135 (1997).
[37] R.F.S. Andrade, D.O. Cajueiro, C.S. Ferreira, Physica A, 295, 323 (2001)
[38] A. Ralston e P. Rabinowitz, A First Course in Numerical Analysis (McGraw-Hill,
Inc. 1978).
[39] C. A. Spindt, I. Brodie, L. Humphrey, E. R. Westerberg, J. Applied Phys, 47,
5248 (1976).
[40] Y. Feng, J. P. Verboncoeur, Phys. Plasmas, 13, 073105 (2006)
105
Apêndice A
Caracterı́sticas e Propriedades de
Alguns Fractais Auto-Similares
Neste apêndice consideraremos alguns fractais que, por sua importância histórica ou riqueza de caracterı́sticas, constituem exemplos ilustrativos “clássicos” de
propriedades de fractais.
A.1
Conjunto de Cantor
A construçcão geométrica do conjunto de Cantor recebe, por vezes, o nome
de “Poeira de Cantor”. Para sua construção inicia-se com um segmento de reta de
comprimento unitário. Divide-se este segmento em 3 partes iguais, retirando-se o
seu terço médio. Essa é a primeira etapa, ou primeiro nı́vel, da construção. Na
segunda etapa, retira-se o terço médio de cada um dos dois segmentos restantes da
primeira etapa. As porções restantes são novamente divididas e delas são retirados os
terços médios, procedendo-se sucessivamente do mesmo modo. O processo é repetido
fazendo-se o número de etapas, ou nı́veis N , tenderem a um número infinitamente
grande. A figura obtida quando N → ∞ é o conjunto de Cantor. Algumas etapas da
sua construção são mostradas na Figura A.1.
É interessante analisar o que acontece com o número de segmentos, nN , com o
comprimento de cada segmento, cN , bem como com o comprimento total do conjunto,
CtN , em cada geração N de sua construção. Entende-se por comprimento total a
soma dos comprimentos dos segmentos de um conjunto. Na nı́vel inicial, ou seja, para
N = 0, tem-se um segmento de modo que n0 = 1. No nı́vel 1, tem-se 2 segmentos. No
nı́vel 2 são quatro segmentos, enquanto na geração 3 são oito segmentos. Deste modo,
106
Figura A.1: Cinco primeiros nı́veis de construção do Conjunto de Cantor.
pode-se provar, por indução finita, que no N − ésimo nı́vel, o número de segmentos
é 2N , ou seja:
nN = 2 N
(A.1)
Prova 2 Seja nN =2N uma propriedade denotada por PN .
- Para N=0, tem-se n0 =1. De fato, ao iniciar-se a construção, o número de
segmentos que se considera é igual a 1. Logo, PN é verdadeira para N=0.
- Seja k ≥ 0. Supondo Pk verdadeira, isto é, nk =2k , deve-se mostrar que Pk+1
também é verdadeira, ou seja, nk+1 =2k+1 .
Naturalmente, por construção, nk+1 = nk .2, pois a cada nı́vel, o número de
intervalos existentes é multiplicado por 2 em relação ao número de intervalos imediatamente anterior. Usando a hipótese de indução, tem-se:
nk+1 = 2k .21 ⇒ nk+1 = 2k+1
(A.2)
Logo, para qualquer N ≥ 0, tem-se nN =2N (CQD).
No Conjunto de Cantor, isto é, para N →∞, tem-se:
lim 2N = ∞
N →∞
(A.3)
Ou seja, no conjunto de Cantor, o número de segmentos tende ao infinito.
Uma contradição com a afirmação anterior pode ser verificada ao se analisar o
comprimento total do conjunto de Cantor, CtN . Para isso, é necessário primeiramente
analisar o comprimento de cada segmento, cN , que compõe o Conjunto de Cantor
no correspondente nı́vel. No primeiro nı́vel, N = 0, o comprimento do segmento é
107
cN =1; no segundo nı́vel cN = 31 ; no terceiro nı́vel cN = 19 . Então, no N-ésimo nı́vel, o
comprimento de cada segmento é expresso por:
µ ¶N
1
3
Portanto tem-se que, no limite para infinitos nı́veis que:
cN =
(A.4)
µ ¶N
1
=0
(A.5)
3
Desta maneira, o comprimento de cada segmento tende a zero. Por isso,
lim cN = lim
N →∞
N →∞
o resultado do conjunto de Cantor é uma série de pontos “pulverizados”; daı́ a
denominação de “Poeira de Cantor”.
Para se analisar o comprimento total CtN do conjunto de Cantor, basta que se
multiplique o número de segmentos pelo comprimento de cada um deles. Logo:
µ ¶N
2
3
Quando o número de gerações tende a infinito, tem-se:
CtN =
(A.6)
µ ¶N
2
=0
(A.7)
N →∞
N →∞ 3
Portanto, o comprimento do Conjunto de Cantor tende a 0. Observam-se pois
lim CtN = lim
caracterı́sticas do conjunto de Cantor que são paradoxais. Ao mesmo tempo que
o número de segmentos de que é composto o conjunto tende a infinito, conclui-se
que o conjunto de segmentos possui um comprimento total nulo. As construções
matemáticas que encerram tais contradições são comumente chamadas de “Monstros
Matemáticos” ou ainda de “Casos Patológicos”.
A.2
Ilha de Koch
Nesta subseção é discutido um dos processos de formação do que se denomina
de ilha de Koch. No caso parte-se de uma linha fechada, denotada de “ilha”, que tem
como geração inicial a forma de um triângulo equilátero. O processo de construção
se inicia substituindo-se o terço central de cada um dos lados, supostos cada um de
108
comprimento unitário, por outros dois segmentos com comprimentos de 13 , formandose uma estrutura triangular, sem a base que justamente é corresponderia à porção
removida. Obtem-se então uma estrutura com comprimento total de 4 unidades (três
conjuntos de quatro partes cada um, cada parte com comprimento 13 ). O processo
é repetido para cada um dos doze segmentos, e assim sucessivamente, até que, para
infinitos nı́veis, tem-se a estrutura chamada de ilha de Koch (ver Figura A.2).
Figura A.2: Os quatro primeiros nı́veis da ilha de Koch triangular.
Primeiramente procede-se a uma análise de como a área, limitada pelos segmentos que constituem a figura, muda no processo iterativo. No nı́vel inicial, N = 0,
tem-se um triângulo equilátero de lado l, cuja área S0 é dada por:
√
l2 3
S0 =
4
109
(A.8)
Para o segundo nı́vel, ou seja, N = 1, tem-se uma ilha limitada por 3x4
segmentos, de modo que a área da ilha S1 será:
Ã !2 √
l
3
(A.9)
3
4
Para o terceiro nı́vel, ou seja, para N = 2, tem-se uma ilha limitada por 3x4x4
S1 = S0 + 3
segmentos, de modo que a área da ilha S2 será:
Ã !2 √
Ã !2 √
l
3
l
S2 = S0 + 3
+ 12
3
4
9
Logo, para o N-ésimo nı́vel, tem-se:
√
3
4
(A.10)
µ ¶2i
N
3X
1
SN = S0 + 3l
4i−1
4 i=1
3
A equação A.11 pode ser reescrita como:
2
√
(A.11)
µ ¶
N
3X
4 i−1
SN = S0 + l .
(A.12)
12 i=1 9
Fazendo N = ∞, tem-se que o somatório da equação A.12 é uma série geométrica
2
de razão q = 49 , de modo que tem-se:
∞ µ ¶i−1
X
4
9
(A.13)
5
i=1 9
Portanto, tomando-se o limite da área SN para infinitos nı́veis, tem-se:
=
2√ 2
3.l
(A.14)
N →∞
5
Este resultado mostra que a área delimitada por uma linha de comprimento
lim SN =
infinito, de acordo com o processo de construção exposto nesta subseção, é finita.
O fato do comprimento da linha, que é fronteira da ilha de Koch, ser infinito, foi
apresentado explicitamente no Capı́tulo 2 para descrever a construção da curva de
Koch.
110
A.3
Triângulo de Sierpinski
Nesta subseção será discutido o processo de construção do Triângulo de Sier-
pinski, assim como algumas caracterı́sticas geométricas do mesmo, como o cálculo da
área obtida para diferentes nı́veis e o cálculo do comprimento resultante para infinitos
nı́veis no processo de construção. Um processo simples de construção do Triângulo de
Sierpinski se inicia a partir de um triângulo equilátero totalmente preenchido (nı́vel
inicial, N = 0). Posteriormente determinam-se os pontos médios de cada um dos
três segmentos que delimitam o triângulo inicial, de modo que ligando-se estes três
pontos médios, obtém-se quatro triângulos cujos lados correspondem à metade do lado
do triângulo inicial. Ao retirar-se o triângulo central tem-se a segunda configuração
correspondente a N = 1, concluindo-se portanto o processo básico de construção. O
processo é repetido com cada um dos três triângulos restantes (Ver Figura A.3), e
sucessivamente com cada triângulo equilátero formado na seqüência.
Figura A.3: Os cinco primeiros nı́veis de construção do Triângulo de Sierpinski.
Para determinação da área do Triângulo de Sierpinski, considera-se inicialmente
111
um triângulo equilátero de lado l, cuja área S0 é dada por:
√
l2 3
S0 =
(A.15)
4
Em cada passo N , subtrai-se a área de nN triângulos com lados lN . Tem-se
então que n1 = 1, n2 = 3, n3 = 9, ..., nN = 3N −1 . Os lados lN , são obtidos pela
redução dos lados do triângulo original por um fator 12 , ou seja:
µ ¶N
1
l
2
Desse modo, a área S1 , será expressa por:
lN =
(A.16)
Ã !2 √
l
2
S1 = S0 −
3
4
(A.17)
A área S2 , será dada por:
Ã !2 √
Ã !2 √
l
3
l
S2 = S0 −
−3
2
4
4
Portanto, para N → ∞, obtém-se:
√
3
4
(A.18)
µ ¶N
∞
3 X
3
SN = S0 − l
12 N =1 4
2
=0
(A.19)
É interessante discutir o perı́metro de cada um dos triângulos obtidos em cada
nı́vel da construção para se calcular a soma do perı́metro dos 3N triângulos no nı́vel N .
Esta análise permite calcular a soma dos perı́metros dos triângulos da figura como um
todo. Analogamente ao que foi feito para determinação do número de triângulos nN ,
determina-se o perı́metro de cada triângulo para um determinado nı́vel. Para N = 0,
tem-se o perı́metro do triângulo original de lado l, ou seja, T0 =3l. Para N = 1, o lado
de cada triângulo gerado será 2l , o que resulta em um perı́metro 3 2l . Para N = 2, o
lado de cada triângulo gerado será 4l , o que resulta em um perı́metro 3 4l , para cada
triângulo. Portanto, para o nı́vel N , tem-se que o perı́metro de cada triângulo, TN ,
resulta em:
TN =
3l
2N
112
(A.20)
Como o número adicional de triângulos removidos do nı́vel N é 3N , a soma PN
dos perı́metros dos triângulos no nı́vel N é dada por:
PN = 3
µ ¶N
3
l
2
Para o Triângulo de Sierpinski, N → ∞, de modo que:
(A.21)
lim PN = ∞
(A.22)
N →∞
Logo, o perı́metro aumenta indefinidamente à medida que aumentamos o número de nı́veis na construção do Triângulo de Sierpinski. Isto leva à uma conclusão
conflitante: a área total de todos os triângulos tende para zero enquanto que o
perı́metro da estrutura formada, aumenta indefinidamente.
A.4
Esponja de Menger
A construção da Esponja de Menger é baseada no mesmo princı́pio utilizado
para construção do Triângulo de Sierpinski, contudo, o processo iterativo é feito com
um cubo, estendendo-se portanto a uma situação tri-dimensional.
O processo de construção se dá de tal forma que para N = 0, tem-se um cubo
maciço de lado l e com volume V0 =l3 . Para N = 1, o cubo é dividido em 27 cubos
menores e iguais, cada um com uma aresta igual a 3l . Remove-se o cubo central, bem
como os seis cubos situados no meio de cada face (Ver figura A.4). Este processo
é repetido seqüencialmente com os todos os cubos restantes, dividindo cada um em
27 outros com
1
3
da aresta do cubo imediatamente anterior. Similarmente, remove-se
o cubo central e cada cubo na porção central das faces. No segundo nı́vel, ou seja,
N = 1, o volume da esponja, V1 , será dado por:
Ã !3
l
(A.23)
3
No terceiro nı́vel, ou seja, N = 2, cada um dos 20 cubos restantes são divididos
V1 = V0 − 7
em mais 27 iguais, dos quais 7 so retirados, cada um com volume
³ ´3
l
9
A.5). Deste modo, o volume da esponja, V2 , será dado pela expressão:
113
(Ver Figura
Figura A.4: Nı́vel N = 1 da Esponja de Menger .
Ã !3
l
V1 = V0 − 7
3
Ã !3
−7
l
9
20
(A.24)
Figura A.5: Geração N = 2 da Esponja de Menger .
Portanto, para o N -ésimo nı́vel e fazendo N →∞, o volume da esponja será
dado por:
∞ µ ¶3N
X
1
20N −1 = 0
(A.25)
3
N =1
Logo, observa-se que o volume da Esponja de Menger tende a zero, quando o
VN = V0 − 7l
3
número de nı́veis tende a infinito. No entanto, para determinação da área da superfı́cie,
SN , com N →∞, desta estrutura fractal, procede-se de tal forma que para N = 0 temse S0 =6l2 ; para N = 1 tem-se:
S1 = S0 + 6l2
114
µ ¶2
1
3
20
(A.26)
Portanto, para o N -ésimo nı́vel, e fazendo N →∞, a área da superfı́cie associada
à enponja será dada por:
∞ µ ¶2N
X
1
20N = ∞
(A.27)
3
N =1
Conclui-se então, que a Esponja de Menger possui volume nulo, bem como uma
SN = S0 + 6l
2
área infinita para infinitos nı́veis.
115
Apêndice B
Equação KPZ
Para o modelo de deposição balı́stica, pode-se determinar os expoentes caracterı́sticos a partir da construção de uma equação contı́nua que na verdade corresponde
à generalização da equação de Edwards-Wilkinson [EW] [4]. Tal modelo se constitui
de correlações laterais e de um crescimento local normal à interface. Para incluir
o crescimento lateral na equação de crescimento, primeiro adiciona-se uma nova
partı́cula à superfı́cie. Para tal crescimento tem-se então, um aumento δh ao longo
do eixo que está associado ao crescimento, que facilmente pode ser determinado pelo
esquema mostrado na Figura B.1. Portanto:
1
1
δh = [(vδt)2 + (vδt∇h)2 ] 2 = vδt[1 + (∇h)2 ] 2
(B.1)
Figura B.1: O esquema acima mostra a origem do acréscimo do termo não linear à
equação [EW] tratando o crescimento local normal à interface.
Se ∇h ¿ 1, pode-se expandir B.1 de modo que:
116
∂h(x, t)
v
= v + (∇h)2 + ...
∂t
2
(B.2)
A equação B.2 , sugere que um termo não linear da forma (∇h)2 deve estar
presente na equação de crescimento, para refletir a presença do crescimento lateral.
Portanto, adicionando o termo acima à equação de EW, obtém-se a equação de
Kardar-Parisi-Zhang (KPZ). Logo:
∂h(x, t)
λ
= v∆h + (∇h)2 + η(x, t)
∂t
2
(B.3)
Na equação acima, o primeiro termo descreve a relaxação da superfı́cie causada
pela tensão de superfı́cie . O ruı́do η(x, t) , satisfaz < η(x, t) >= 0 e < η(x, t)η(x0 , t0 ) >=
2Dδ d (x−x0 )δ(t−t0 ) , sendo que a primeira reflete as flutuações randômicas no processo
de deposição com número randômico descorrelacionado de média nula e a segunda
identifica que o ruı́do não apresenta correlações no espaço e no tempo desde que a
média do produto dos ruı́dos seja nula, exceto para o caso especial em que x = x0 e
t = t0 , sendo d a dimensão do substrato e D uma constante.
Outra caracterı́stica importante da equação KPZ, é que a mesma apresenta
simetrias associadas à invariância de translação no tempo, ou seja, não depende de
onde se define a origem do tempo sendo portanto invariante à transformação t → t+δt
. Uma segunda simetria presente, seria ao longo da direção de crescimento, sendo
portanto independente de onde se define h = 0 , ou seja, é invariante à transformação
h → h + δh. Porém, a simetria da altura de interface é quebrada, uma vez que as
mesmas estão intimamente ligadas à natureza de equilı́brio da interface formada, ou
seja, no processo de deposição balı́stica a propriedade de crescimento lateral existe,
tendo como conseqüência a não invariância da equação KPZ diante da transformação
h → −h. É importante ressaltar que o termo não linear afeta os expoentes de escala
associados à interface auto-afim.
Para tal modelo, as transformações de escala que tornam a interface estatisticamente indistinguı́vel são:
117
x → x0 ≡ bx
h → h0 ≡ bα h
t → t0 ≡ bz t
Substituindo-se as transformações acima na equação B.3 e usando-se a propriedade da função delta δ d (ax) =
bα−z
1
δ(x)
ad
, tem-se:
1
∂h
λ
= vbα−2 ∆h + b2α−2 (∇h)2 + b− 2 (d+z) η(x, t)
∂t
2
(B.4)
Se compararmos bα−2 ∇2 h com b2α−2 (∇h)2 , nota-se que no limite b → ∞, o
termo não linear domina o termo que envolve a tensão de superfı́cie, com α > 0 .
Multiplicando-se ambos os lados da equação B.4 por bz−α , então:
1
∂h
λ
= vbz−2 ∆h + bα+z−2 (∇h)2 + b− 2 (d−z)−α
∂t
2
(B.5)
Um detalhe importante da equação acima, na determinação dos expoentes
caracterı́sticos é que os termos v, λ e D na equação de crescimento não se renormalizam
independentemente sendo, portanto, acoplados. Então, para um comportamento de
invariância de escala, não se pode simplesmente assumir que os expoentes de b são
nulos, pois os coeficientes v, λ e D podem também mudar após a transformação de
escala.
118
Apêndice C
Método de Liebmann
A solução analı́tica de uma equação diferencial nem sempre é possı́vel. No caso
de equações diferenciais parciais de segunda ordem, esta dificuldade geralmente tem
origem no fato de que a solução buscada, u = u(x, y) (no caso de duas dimensões),
deve ser válida em uma região do espaço com fronteira definida, sobre a qual uma
condição de contorno será imposta. A topologia desta fronteira, se irregular, constitui
dificuldade à possibilidade de uma abordagem analı́tica. Numa situação deste tipo,
resulta inevitável recorrer a uma metodologia numérica.
O emprego de técnicas numéricas para a solução de equações diferenciais parciais, requer que a região de integração, em lugar de ser tratada como um contı́nuo,
constitua um conjunto discreto e finito de pontos nos quais a variável de interesse é
calculada. Quanto maior for o número de pontos em que a região é “dividida”, mais
próxima da solução exata será a solução aproximada. O conjunto de pontos discretos
constitui assim uma malha. Uma vez que o domı́nio é tratado discretamente, obtemse um conjunto de equações, onde o valor da variável num ponto é escrito em função
dos valores da mesma variável em outros pontos da malha. Este procedimento resulta
em um sistema de equações algébricas, geralmente lineares, que representa a equação
diferencial parcial. O procedimento é iterativo de modo que a solução convirja.
Os métodos iterativos partem de uma aproximação inicial, digamos u(0) , e
constróem uma sucessão de soluções, u(1) , u(2) , ..., u(k) que, espera-se, seja convergente
para a solução exata u do problema. Seja Au = B um sistema linear, onde A é uma
matriz m x n, u um vetor coluna de dimensão n e B uma matriz linha de dimensão m.
Converte-se sistema num equivalente do tipo u = M u+c, onde M é uma matriz m x n
119
e c é um vetor coluna de dimensão n. Assim, o método iterativo, do tipo estacionário,
pode ser escrito da seguinte forma:
u(k+1) = M u(k) + c
(C.1)
em que M é designada matriz iteração.
Dentre as equações diferenciais parciais elı́pticas, a equação de Laplace pode ser
resolvida usando a técnica de diferenças finitas, transformando-se em uma equação
algébrica de diferenças. As condições de fronteira exigidas para a solução podem
ser especificadas atribuindo-se o valor da variável u, cuja solução se procura, na
fronteira da região de integração. Este tipo de condição é chamado de condição de
fronteira de Dirichlet. Alternativamente, as primeiras derivadas da função u podem
ser especificadas na fronteira, o que neste caso se chama de condição de Neumann.
A equação diferencial parcial é resolvida dividindo-se o domı́nio de integração
em uma malha (“grid”) retangular e escrevendo a equação de diferanças para cada
um dos pontos que compõem a malha, denotados por nodos. Isto leva a um conjunto
de equações algébricas lineares que podem ser resolvidas simultaneamente.
Na argumentação que se segue, assume-se que a função que se deseja determinar
seja o potencial elétrico. Os valores do potencial em cada ponto são inicialmente
desconhecidos e determinados iterativamente utilizando-se um critério de convergência.
Assim, repete-se o processo de cálculo até que o vetor u(k) , em cada ponto da malha,
possa ser considerado como suficientemente próximo de u(k−1) , ou seja, define-se uma
tolerância ε de modo que:
°
°
° (k)
°
°u − u(k−1) ° ≤ ε
(C.2)
Seja h = ∆x o espaçamento do “grid” na direção x, e a função u(x) como sendo
contı́nua até a quarta derivada. Utilizando-se uma expansão em série de Taylor, para
o caso bi-dimensional, tem-se:
120
1
1
1
u(xn +h) = u(xn )+u0 (xn ).h+ u00 (xn ).h2 + u000 (xn ).h3 + u(iv) (ξ1 ).h4 , xn < ξ1 < xn +h
2!
3!
4!
(C.3)
e portanto,
1
1
1
u(xn −h) = u(xn )−u0 (xn ).h+ u00 (xn ).h2 − u000 (xn ).h3 + u(iv) (ξ2 ).h4 , xn −h < ξ2 < xn
2!
3!
4!
(C.4)
Somando-se as equações C.3 e C.4 tem-se:
u(xn + h) − 2u(xn ) + u(xn − h)
uiv (ξ) 2
00
=
u
(x
)
+
.h , xn − h < ξ < xn + h (C.5)
n
h2
12
Podemos simplificar a notação da equação C.5 indicando, como ı́ndice de u, n
que se refere à posição do ponto no grid e por ± indicando se u é calculado em x ou
num ponto em relação a x deslocado de ±h. Já o termo de ordem O(h2 ) significa que
o erro é proporcional a h2 quando h → 0. Portanto:
un+1 − 2un + un−1
' u00n + O(h2 )
(C.6)
h2
De modo similar, a primeira derivada pode ser aproximada como dependendo
da diferença entre os valores das funções nos pontos sucessivos do “grid”, tendo-se
portanto:
un+1 − un−1
= u0n + O(h2 )
(C.7)
2h
Com as aproximações acima, no ponto (xi , yi ) o Laplaciano assume a forma:
ui+1,j − 2ui,j + ui−1,j ui,j+1 − 2ui,j + ui,j−1
+
(C.8)
(∆x)2
(∆y)2
Como, no caso, desejamos resolver a equação de Laplace, ∆u(xi , yi ) = 0.
∆u(xi , yi ) '
Considerando-se um “grid” regular, com h = ∆x = ∆y = 1, que substituindo na
eq. C.8 leva a que o potencial em cada ponto seja dado por:
121
1
ui,j = {ui+1,j + ui−1,j + ui,j+1 + ui,j−1 }
(C.9)
4
Na equação C.9 temos os valores do potencial em cinco pontos distintos: à
diteita, à esquerda, acima e abaixo de um ponto “central” ui,j (Ver Figura C.1). Para
a determinação do potencial em todos os pontos, temos que distinguir pontos nos
quais o valor do potencial varia ao longo das iterações e pontos para os quais o valor
do potencial é mantido constante ao longo de todo o processo iterativo. Estes últimos
correspondem aos pontos que definem a fronteira, portanto estamos utilizando uma
condição de contorno tipo Dirichlet. Valores iniciais são arbitrariamente atribuı́dos
a cada um dos pontos do “grid”, que não os de fronteira, obedecendo à condição
de que estes valores sejam compatı́veis com os potenciais previamente atribuı́dos aos
contornos. Assim, na execução computacional, a cada ponto do domı́nio de integração,
associa-se uma propriedade, que pode ser chamada de função estado, atribuı́do o valor
1 para pontos onde o potencial é mantido inalterado e 0 para pontos onde o potencial
é recalculado pelo procedimento iterativo. Mantendo-se constante o potencial da
fronteira, o processo de cálculo iterativo do potencial em cada ponto é realizado
utilizando-se a eq. C.9, até que a diferença entre os potenciais calculados, num mesmo
ponto e para todos os pontos, em sucessivas iterações, se torne menor do que um valor
adotado previamente, escolhido como um critério de convergência. Há a possibilidade
de que a fronteira da região onde desejamos determinar o potencial, corresponda
a dois ou mais conjuntos de pontos onde o potencial, distinto para cada conjunto,
seja mantido constante. É possı́vel, por exemplo, que a região tenha como fronteira
duas superfı́cies mantidas a uma diferença de potencial. Neste caso, os extremos
do “grid”, numa direção, definem o que poderı́amos chamar de “laterais” da região
de integração. Com relação às laterais pode-se adotar, por exemplo, condições de
periodicidade evitando portanto a atribuição de potenciais definidos.
No caso tri-dimensional utiliza-se, para se escrever o laplaciano, três indices
que indicam as posições espaciais dos pontos. Para um “grid” genérico, tem-se:
122
Figura C.1: Região onde se deseja determinar o potencial elétrico, tratada como um
conjunto discreto de pontos.
∆u =
ui+1,j,k − 2ui,j,k + ui−1,j,k ui,j+1,k − 2ui,j,k + ui,j−1,k ui,j,k+1 − 2ui,j,k + ui,j,k−1
+
+
(∆x)2
(∆y)2
(∆z)2
(C.10)
Sendo ∆u(xi , yi , zi ) = 0, o potencial em cada ponto do “grid” será dado por:
ui,j,k
(ui−1,j,k + ui+1,j,k ).L2Ψ .L2ξ (ui,j−1,k + ui,j+1,k )L2θ .L2ξ (ui,j,k−1 + ui,j,k+1 )L2Ψ .L2θ
=
+
+
2(L2Ψ .L2ξ + 1 + L2Ψ )
2(L2θ .L2ξ + 1 + L2ξ )
2(L2Ψ .L2θ + 1 + L2θ )
(C.11)
onde LΨ =
∆y
,
∆z
Lξ =
∆z
∆x
e Lθ =
∆x
∆y
e onde ∆x, ∆y e ∆z correspondem aos
espaçamentos entre pontos vizinhos nas direções x, y e z, respectivamente. É possı́vel
notar que a regularidade do grid, ou seja, no caso em que ∆x = ∆y = ∆z, reduz
a eq.C.11 a uma média aritmética do potencial nos pontos que correspondem aos
seis primeiros vizinhos de um ponto central. O conjunto de equações necessário para
este caso é mais extenso mas, em princı́pio, os métodos são idênticos aos expostos na
solução para o caso bi-dimensional.
Em problemas tri-dimensionais, torna-se fácil exceder o espaço de memória
computacional disponı́vel. Por exemplo, se o volume considerado tem 100 pontos em
cada direção, um total de (100)3 pontos estarão envolvidos. Isto pode requerer a
123
representação dos coeficientes de um milhão de equações numa matriz quadrada.
124
Apêndice D
Algoritmo de Verlet
O algoritmo de Verlet é provavelmente o método de integração mais utilizado
em dinâmica molecular. O mesmo utiliza a expansão em série de Taylor, levando-se
em conta os termos de terceira ordem da posição da partı́cula em função do tempo.
Deste modo, expandindo-se o vetor posição em um instante t + ∆t, tem-se que:
1
1
(D.1)
~rj (t + ∆t) = r~j (t) + ~vj (t)∆t + ~aj (t)(∆t)2 + ~bj (t)(∆t)3 + O(∆t)4
2
6
Na equação D.1, ~aj (t) = d~vj (t) é a aceleração e o vetor ~bj (t) corresponde à
dt
derivada da aceleração. Aplicando-se, à referida equação, o sentido de evolução
inverso, tem-se:
1
r~j (t − ∆t) = r~j (t) − ~vj (t)∆t + ~aj (t)(∆t)2 −
2
1~
bj (t)(∆t)3 + O(∆t)4
6
(D.2)
Somando-se as equações D.1 e D.2, chega-se à equação fundamental do algoritmo,
ou seja :
r~j (t + ∆t) + r~j (t − ∆t) = 2~
rj (t) + ~aj (t)(∆t)2 + O(∆t)4
(D.3)
Uma vez que as coordenadas cartesianas de cada partı́cula são tratadas como
coordenadas generalizadas e o potencial de interação depende exclusivamente das
posições, ou seja, U = U (~
r1 , ..., r~N ), tem-se pela equação de Newton que:
mj r~¨j = F~j
(D.4)
Como tem-se que as forças são derivadas do potencial, a força total exercida
125
sobre a j-ésima partı́cula, F~j , pode ser expressa por:
~ r U (~
F~j = −∇
r1 , ..., r~N )
j
(D.5)
Substituindo-se o valor da aceleração da equação D.4 em D.3, tem-se que:
r~j (t + ∆t) + r~j (t − ∆t) = 2~
rj (t) +
Fj (~
r1 (t),~..., r~N )(t)
(∆t)2 + O(∆t)4
mj
(D.6)
Observa-se portanto que o algoritmo de Verlet possui simplicidade de implementação, sendo estável e com boa precisão uma vez que seu erro de truncamento
é da ordem de (∆t)4 . Contudo, o mesmo não calcula automaticamente a velocidade
uma vez que esta última não é necessária para a determinação das trajetórias. Para se
determinar grandezas, como por exemplo a energia cinética, K, que pode ser utilizada
para se comprovar a conservação da energia mecânica, E = K + U , determina-se a
velocidade da partı́cula subtraindo-se as equações D.1 e D.2 de modo que:
r~j (t + ∆t) − r~j (t − ∆t)
+ O(∆t)2
(D.7)
2∆t
Vale ressaltar que, devido à sua simplicidade, o algoritmo de Verlet requer um
~vj (t) =
reduzido recurso computacional, fato este que não chega a apresentar uma grande
vantagem.
Tal algoritmo requer também que o intervalo de tempo utilizado na
simulação seja constante. Porém, a variação do ∆t no decorrer da simulação pode levar
a uma vantajosa economia no tempo de simulação. Pode-se, por exemplo, recorrer a
situações onde as forças entre as partı́culas sejam grandes implicando em acelerações
altas. Como conseqüência, os intervalos de tempo teriam de ser menores. Quando as
forças são de baixa intensidade pode-se empregar intervalos de tempo maiores. Isto
sugere o ajuste do intervalo de tempo, ∆t, conforme o comportamento do sistema, o
que possibilita otimizar o tempo de processamento consideravelmente.
126