METODOLOGIA PARA ESTUDO COMPUTACIONAL DO DESEMPENHO HARMÔNICO DE SISTEMAS ELÉTRICOS DE GRANDE PORTE LEANDRO R. ARAUJO 1, SERGIO L. VARRICCHIO2, SERGIO GOMES JR2, NELSON MARTINS2, CRISTIANO O. COSTA 1, FRANKLIN C. VÉLIZ 1 1 – FCT-JF – Fundação Centro Tecnológico, Juiz de Fora, MG, Brasil 2 – CEPEL, Centro de Pesquisa em Energia Elétrica Caixa Postal 68007, Rio de Janeiro, RJ, Brasil E-mail:
[email protected],
[email protected],
[email protected],
[email protected],
[email protected],
[email protected] Abstract This paper presents new methodologies and tools for the power systems harmonics studies in large scale systems. The proposed tools follow the paradigm of the Object Oriented Modeling; and make possible the connection of several subsystems. Examples showing the advantages of using the proposed tools for the harmonic studies are presented. The methodologies were implemented in HarmZs, a production-grade program for harmonic studies in electric power systems. Keywords Harmonics, Electrical Network Modeling, Object Oriented Modeling, Large Scale Systems, Current Injection Methods Resumo Este artigo apresenta novas metodologias e ferramentas para o estudo computacional do desempenho harmônico para sistemas elétricos de grande porte. As ferramentas propostas seguem o paradigma da Modelagem Orientada a Objetos, e possibilitam a conexão de diversas sub-redes. Também são apresentados exemplos que mostram as vantagens da utilização das ferramentas propostas para os estudos harmônicos. As metodologias foram implementadas no HarmZs, um programa de nível de produção para estudos harmônicos em sistemas elétricos de potência. Palavras-chave Harmônicos, Modelagem de Redes Elétricas, Modelagem Orientada a Objetos, Sistemas de Grande Porte, Método de Injeção de Corrente.
1
Introdução
O programa HarmZs [1] é uma ferramenta computacional baseada em tecnologias mais recentes, tanto na modelagem de redes elétricas quanto nas técnicas de programação computacional. As tecnologias de modelagem de redes elétricas são denominadas matriz Y(s) [2]-[5] e Sistemas Descritores [5]-[8]. Ambas permitem a análise harmônica em todo o plano complexo s = σ + jω, sendo a possibilidade da análise modal a principal vantagem deste domínio estendido. As versões do programa HarmZs são desenvolvidas utilizando Microsoft Visual C++ 6.0 [9]-[10] e Microsoft Foundation Classes (MFC) [10], pois a MFC possui muitas classes implementadas que facilitam o desenvolvimento de aplicativos compatíveis com o padrão Windows. O HarmZs também segue o paradigma de desenvolvimento de programas orientado a objetos [11]-[14] e a tecnologia de componentes por meio do Active X e do Component Object Model (COM) [15]-[16]. A programação gráfica permite uma interface amigável com o usuário. Os resultados das simulações apresentados em gráficos e/ou tabelas tornam suas análises mais fáceis e rápidas. Alguns conceitos básicos da análise harmônica convencional [17] são descritos neste artigo, pois são necessários para a utilização do programa HarmZs, são revisados. Algumas implementações computacionais que facilitam a aplicação do programa a sis-
temas de grande porte são apresentadas. Dentre elas destaca-se uma nova metodologia para a configuração de redes elétricas a partir de sub-redes. Como exemplo de aplicação desta nova metodologia, um estudo de distorção harmônica em um sistema industrial é apresentado considerando duas situações: 1. O sistema de transmissão (sistema brasileiro) ao qual a indústria em análise está conectada é representado simplesmente por sua impedância de Thévenin. 2. O sistema de transmissão é representado integralmente com todos os seus componentes. 2
Metodologias Convencionais
As metodologias convencionais para análises de comportamento harmônico de sistemas de potência podem ser divididas em dois grupos, de acordo com a modelagem adotada para a rede elétrica: 1. Domínio do tempo 2. Domínio da freqüência Como exemplo de ferramentas computacionais que pertencem ao grupo 1, podem-se citar os programas EMTP, ATP e PSCAD que são baseados na integração trapezoidal das equações dos elementos de rede escritas no domínio do tempo. Como exemplo de metodologia que pertence ao grupo 2 pode-se citar o método de injeção de correntes.
1
A equivalência entre os dois métodos pode ser verificada com o auxílio do circuito apresentado na Figura 1. Este circuito contém um reator controlado a tiristores (RCT), os quais fazem com que o reator conduza corrente em determinados intervalos de tempo. O comprimento destes intervalos é determinado pelo ângulo de disparo (α) dos tiristores, medido a partir do cruzamento da tensão da barra 2 por zero. R12
1
L12
onde a0 / 2 é igual ao valor médio da função f(t), ω0 a freqüência fundamental e Dn e φn representam o módulo (valor eficaz) e o ângulo do n-énesio harmônico. Suas expressões analíticas são dadas por Dn =
b a n2 + bn2 e φ n = − atan n 2 an
sendo an =
2
bn = v2
vf
C2 iL2
T
2 T
∫ f (t ) cos (n ω t ) dt
2 T
T
(3)
0
0
∫ f (t ) sen (n ω t ) dt
(4)
0
0
Os coeficientes Dn e os ângulos φn dos regimes permanentes da tensão na barra 2 (Figura 2) e da corrente no elemento L2 (Figura 3) estão mostrados na Tabela 2.
L2
Figura 1: Circuito com RCT. 80
2.1 Análise no Domínio do Tempo
60
Tabela 1: Dados do Circuito da Figura 1
α 1200
R12 (Ω) L12 (H) C2 (µF) L2 (H) 291.75 3.8661 0.35014 12.057
40
Corrente (A)
O circuito mostrado na Figura 1 foi modelado no programa PSCAD com os valores dos parâmetros listados na Tabela 1.
20 0 -20 -40
∆t (µs) 10
-60 -80 0
0.04
0.08
A fonte de tensão é descrita pela expressão (1). vf =
2
500 sen (ωt ) kV
3
n
Tensão (kV)
400 200 0 -200 -400 -600 0
0.04
0.08
0.12
0.16
0.2
Tempo (s)
Figura 2: Tensão v2
Devido à natureza periódica dos regimes permanentes destas curvas, os mesmos podem ser decompostos em termos de séries de Fourier, cuja teoria é brevemente revista a seguir. Seja f(t) uma função periódica do tempo t de período Τ. A expansão em série de Fourier desta função é dada por f (t ) =
a0 + 2 2
∞
∑ D cos (n ω t + φ ) n
n =1
0
n
(2)
0.2
Tabela 2: Módulo e ângulo dos fasores hamônicos de tensão e corrente.
Na Figura 2 é apresentada a curva de tensão na barra 2 e na Figura 3 a curva de corrente no indutor L2. 600
0.16
Figura 3: Corrente iL2.
(1)
800
0.12
Tempo (s)
1 3 5 7 9 11 13 15
v2 Dn (kV) 296.85 74.427 2.7776 0.83830 0.82257 0.26749 0.077750 0.16647
iL2 φn (graus) -91.116 -86.058 78.787 86.608 -99.949 73.733 79.845 -106.33
Dn (A) 33.508 12.570 1.4607 0.70166 0.90650 0.38017 0.13830 0.31415
φn (graus) 179.36 178.74 -11.520 -3.1427 169.64 -15.568 -8.3262 162.81
2.2 Análise no Domínio da Freqüência As componentes harmônicas de tensão e corrente, obtidas anteriormente resolvendo o circuito mostrado na Figura 1 no domínio do tempo, podem, também, serem determinadas resolvendo o circuito mostrado na Figura 4 no domínio da freqüência. Este circuito é construído considerando os componentes lineares do circuito original, mostrado na Figura 1. O componente não-linear é substituído por uma fonte de corrente com o mesmo conteúdo harmônico da corrente do elemento não-linear (iL2), com a devida inversão de sinal. Uma vez que a fonte de
2
i12 R12
L12
2
v2
30
Distorção de Tensão (%)
tensão é não-nula apenas para 60 Hz, a mesma é representada por um curto-circuito para as demais freqüências harmônicas.
PSCAD
20 15 10 5 0
i2 = -iL2
C2
HarmZs
25
180
300
420
540
660
780
900
Freqüência (Hz)
iC2
Figura 5. Distorção de Tensão na Barra 2. Figura 4. Circuito Harmônico Equivalente.
120
z 2 (s ) = 1 y 2 (s )
(5)
onde: y2 (s ) =
1
z12 (s )
+
1
zC 2 (s )
30 0 -30 -60 -90
(6)
-120 180
As tensões e as correntes do circuito podem, então, ser determinadas para qualquer freqüência de interesse, ou seja: v2 (s ) v (s ) e i12 (s ) = 2 z12 (s ) zC 2 (s )
A distorção individual do harmônico de ordem n na barra k é definida como um valor percentual da tensão base desta barra, ou seja: =
vkn vk (base )
× 100 % , n = 1, 2, K , nmax
300
420
540
660
780
900
Freqüência (Hz)
Deve-se observar que para redes de grande porte, não é possível a determinação analítica das diversas impedâncias nodais. Neste caso, o programa HarmZs, monta a matriz de admitâncias nodais e, por meio de fatoração LDU e solução de vetores esparsos (técnicas de esparsidade), obtém os elementos de interesse da matriz de impedâncias nodais.
d kn
PSCAD
60
z12 (s ) = R12 + s L12 e zC2 (s ) = 1 (s C2 )
v2 (s ) = z 2 (s ) i2 (s ) , iC 2 (s ) =
HarmZs
90
Ângulo da Tensão (%)
A impedância vista dos terminais na fonte i2, ou seja, vista do nó 2 (impedância nodal) é dada por:
(7)
onde nmax denota o número total de harmônicos considerados, ω0 a freqüência fundamental em rad/seg,
v kn o módulo da tensão harmônica de ordem n na barra k e v k ( base ) a tensão base da barra k.
Figura 6. Ângulo de Tensão na Barra 2
Na prática, os valores das componentes harmônicas da corrente injetada também podem ser obtidos utilizando-se um dos seguintes procedimentos: 1.
Considerando a tensão da barra onde o equipamento está conectado como senoidal, ou seja, não é considerada a interação harmônica entre corrente e tensão. 2. Utilizando um programa de fluxo harmônico. 3. Por meio de medições. Para que os dois métodos forneçam os mesmos resultados, a fonte de corrente do método no domínio da freqüência deve possuir o mesmo conteúdo harmônico da corrente do elemento não linear do método do domínio do tempo, pois neste caso a interação harmônica entre tensão e corrente está devidamente considerada. No entanto quando esta interação não é considerada (como, por exemplo, no procedimento 1 descrito acima), os resultados obtidos pelo método no domínio da freqüência serão aproximações dos obtidos pelo método no domínio do tempo. 3
Sistemas de Grande Porte
Para o sistema em questão, tem-se: vk (base ) = 500
3 = 288.68 kV
(8)
Os resultados obtidos por este método são comparados com os obtidos pela integração numérica das equações dos elementos de rede escritas no domínio do tempo, como pode ser observado nas Figuras 5 e 6. Desta observação, concluí-se que os resultados obtidos pelas duas metodologias são praticamente os mesmos.
Algumas das implementações realizadas que possibilitam a utilização do programa para sistemas de grande porte são descritas nesta seção. 3.1 Leitura e Tratamento de Dados Para facilitar a utilização de um programa, é interessante que o mesmo possa retirar informações de bases de dados já existentes e as compatibilize com
3
suas aplicações específicas. Para isto, criou-se no HarmZs um módulo orientado a objetos para aproveitar dados de componentes de redes elétricas disponíveis em arquivos de dados originalmente escritos para os programas ANAREDE [18] e ANATEM [19] do CEPEL. Assim, diversos dados referentes a elementos "shunts", linhas de transmissão, transformadores de 2 e 3 enrolamentos e cargas são aproveitados de arquivos históricos do ANAREDE. Deve-se observar que este aproveitamento de dados não se resume apenas a sua simples leitura, uma vez que diversos dados devem ser tratados para uma correta modelagem da rede para estudos de comportamento harmônico. Por isto, este módulo de importação de dados possui vários métodos com heurísticas como, por exemplo, para identificar possíveis transformadores de três enrolamentos e filtros harmônicos, para calcular parâmetros nominais de linhas de transmissão, para modelagem de cargas em configurações série, paralela ou mista, para tratamento de elementos não-lineares, etc. As informações para a modelagem das máquinas síncronas do sistema são obtidas por leitura de arquivos de dados do programa ANATEM. Na Figura 7 estão mostradas as respostas em freqüência do módulo da impedância de uma barra de 230 kV do sistema elétrico brasileiro, com e sem a modelagem das máquinas síncronas. Como pode ser observado, as diferenças entre as curvas são significativas. Na Figura 8 está apresentado o diálogo principal da interface gráfica do módulo de aproveitamento de dados. Observa-se que neste diálogo existe a opção de consideração de possíveis filtros harmônicos, pois, em geral, no ANAREDE estes filtros são modelados apenas como reatâncias capacitivas na freqüência fundamental. Este tratamento é correto para estudos de fluxo de potência, mas completamente inadequado para estudos de comportamento harmônico. Ao se pressionar o botão Possíveis Filtros o programa apresentará na tela o diálogo mostrado na Figura 9. Neste diálogo o usuário poderá escolher como será feito o tratamento dos possíveis filtros harmônicos. Destacase a terceira opção onde os dados de possíveis filtros são desconsiderados do arquivo principal, mas gravados em um outro arquivo (sub-rede), que pode ser modificado pelo usuário e incorporado futuramente aos dados de rede. 0,10
Considerando as máquinas Desconsiderando as máquinas
|Z| (ohms)
0,07
0,03
0,00 0
500
Freqüência (Hz)
1000
Figura 7. Resposta em Freqüência.
1500
3.2 Sub-Redes Uma grande dificuldade atual em estudos de engenharia envolvendo sistemas elétricos de grande porte, é a impossibilidade, apresentada pela maioria dos programas computacionais, de se montar uma rede elétrica por meio da conexão de sub-redes. Muitas vezes, nestes estudos, de forma a se evitar um grande trabalho manual de construção e compatibilização de arquivos de dados, uma das redes é representada apenas por seu equivalente de Thèvenin.
Figura 8. Conversão de Dados.
Figura 9. Tratamento de Filtros.
Considerando, por exemplo, um estudo de distorções harmônicas em um sistema industrial, a representação da rede de transmissão, na qual a industria está conectada, por sua impedância de Thèvenin na freqüência fundamental pode levar a resultados errados, uma vez que freqüências muito mais elevadas são consideradas neste tipo de estudo. Neste caso, o ideal seria a completa representação da rede de transmissão e do sistema industrial, ou, alternativamente, a representação da rede de transmissão por seu equivalente dinâmico. Com o objetivo de facilitar a completa representação das redes, foi criada uma metodologia que possibilita a conexão de várias sub-redes, praticamente sem esforço adicional. Ou seja, as sub-redes são independentes, podendo ser utilizadas nos estudos separadamente ou conectadas. Para que as sub-redes sejam devidamente conectadas, o usuário necessita indicar apenas as barras equivalentes em cada subrede, sem precisar se preocupar se barras de uma sub-rede possuem a mesma numeração de barras de outras sub-redes. Analogamente, o processo pode ser efetuado para a exclusão de sub-redes. A metodolo-
4
gia não apresenta número máximo de sub-redes a serem conectadas. Uma facilidade criada a partir desta idéia de conexão é a montagem de sub-redes de equipamentos como, por exemplo, de filtros harmônicos, que podem ser conectados facilmente às redes, para serem utilizados apenas em determinados estudos ou casos. Assim, os mesmos são incorporados às redes sem modificar a base de dados original. 3.3 Algoritmo de Conexão de Sub-Redes O modelo de objetos implementado no HarmZs para a configuração de redes elétricas é similar ao modelo apresentado em [14]. Considere as duas sub-redes apresentadas na Figura 10. As mesmas são armazenadas internamente conforme apresentado na Figura 11. Observa-se na Figura 11 que os objetos instanciados apontam-se mutuamente. 1 3
1
2
Barra1 Sub 1
Barra3 Sub 1
Barra2 Sub 1
Linha 1-3
Barra1 Sub 2
Linha 2-3
Barra2 Sub 2
Linha 1-2
Figura 13. Modificações Internas.
3.4 Diagrama de Classes O diagrama de classes que possibilita as operações descritas neste artigo com as redes elétricas é apresentado na Figura 14. Observa-se a existência da classe CRede (CConfigurador), na qual seu objetivo principal é fazer a extração de dados de vários programas do CEPEL através do auxílio da classe CConvCEPEL, ler arquivos de redes no formato HarmZs e reconfigurar as redes com a ajuda da classe CSubRede. O programa HarmZs conecta-se diretamente a classe CRede para a utilização de dados. HarmZs
CConfigurador
2
CConvCEPEL
CSubRede
CIOAnarede
CIOAnatem
CIOFlupot
RET CAS (DLL)
Sub - Rede 1
Sub - Rede 2
Figura 14. Diagrama de Classes.
Figura 10. Sistemas Exemplos. Barra1 Sub 1
Barra2 Sub 1
Linha 1-3
Barra3 Sub 1
Barra1 Sub 2
Linha 2-3
Barra2 Sub 2
Linha 1-2
Figura 11. Armazenamento Interno.
Caso seja feita a equivalência entre a barra 3 da sub-rede 1 com a barra 1 da sub-rede 2, tem-se a rede apresentada na Figura 12. 1 3
4
Estudo de Casos
O sistema industrial a ser estudado é apresentado na Figura 15. O sistema possui 5 barras, filtros para 3º, 5º, 7º e 11º harmônicos. Estes filtros, neste exemplo, foram inseridos como sub-redes. Na barra 5 encontra-se uma fonte de corrente harmônica nas freqüências de ordem 3, 5, 7, 9, 11 e 13. Deve-se observar que neste exemplo existe apenas uma fonte de corrente harmônica. No entanto, diversas fontes de corrente harmônica podem ser consideras [2]. O ponto de conexão desta indústria com o sistema de transmissão é a barra 1.
2
T ransmissão 1
2
2 Figura 12. Barra Equivalentada.
As modificações feitas pelo configurador de redes do HarmZs são apresentadas na Figura 13. Como pode ser observado, criou-se uma nova ligação (tracejada em negrito) da linha 1-2 da sub-rede 2 com a barra 3 da sub-rede 1. A barra 1 da sub-rede 2 (cinza) é armazenada em outra parte e continua apontando para a linha 1-2 da sub-rede 2. O tempo necessário para a completa configuração de redes com aproximadamente 5000 barras é de centésimos de segundos.
3
4 ( AT/MT)
5 Ih1
3
5
7
11
Filtros
Figura 15. Sistema Industrial.
Primeiramente, o sistema de transmissão foi representado por uma impedância de Thèvenin e calculou-se respostas em freqüência do módulo da impe-
5
dância transferência entre as barras 1 e 5 e a distorção harmônica na barra de conexão (barra 1). Depois, os cálculos foram refeitos utilizando-se os dados do ANAREDE e do ANATEM para montar a rede básica e conectou-se o sistema industrial e os respectivos filtros, conforme apresentado na seção 3. Os resultados são apresentados nas Figuras 16 e 17.
redes e de importação de dados dos programas ANAREDE e ANATEM. Estas rotinas computacionais foram implementadas utilizando os conceitos Modelagem Orientada a Objetos. Referências Bibliográficas [1]
0,010
Equivalente (5,1) Completa (5,1)
[2]
|Z| (ohms)
0,007
0,003
[3] 0,000 0
500
1000
1500
Freqüência (Hz)
[4] Figura 16. Impedância de Transferência 1 - 5. [5]
0,45
Equivalente
0,40
Completo Distorção da Tensão (%)
0,35
[6]
0,30 0,25 0,20 0,15 0,10
[7]
0,05 0,00 180
300
420
540
660
780
Freqüências (Hz)
Figura 17. Distorção Harmônica na Barra 1.
Como pode-se observar, as diferenças entre as curvas são significativas, mostrando que não se deve utilizar equivalentes estáticos para representar sistemas em estudos harmônicos.
[8]
[9] [10]
Conclusões Foi inicialmente apresentado um exemplo tutorial de distorção harmônica em um sistema de duas barras, com uma componente não-linear, utilizando-se metodologias de análise no tempo e também na freqüência. Neste exemplo pode-se concluir que os resultados obtidos pelas duas metodologias são praticamente os mesmos. Outro exemplo mostrando a necessidade de se modelar corretamente as máquinas síncronas em um estudo harmônico foi apresentado, verificando-se que desprezar as máquinas acarreta em diferenças significativas. Em outro exemplo foi mostrada a perda de qualidade dos resultados quando o sistema de transmissão, ao qual uma indústria estava conectada, foi representado simplesmente por sua impedância de Thèvenin. Este artigo também apresenta brevemente a metodologia utilizada pelo programa HarmZs para a configuração de redes elétricas. Em particular foram descritos os mecanismos de gerenciamento de sub-
[11] [12]
[13]
[14]
[15] [16] [17] [18] [19]
Varricchio, S. L.; Gomes Jr., S.; Martins, N.; Araujo, L. R. e Costa, C. O. “Manual do Usuário do Programa HarmZs V. 1.2”, CEPEL Relatório Técnico DPP/POL 1385/02, Janeiro 2003. Varricchio, S. L.; Gomes Jr., S.; Martins, N.; “Modal Analysis of Industrial System Harmonics Using the sDomain Approach”, IEEE Transmission and Distribution Latin America Conference, São Paulo, SP, Brasil, 18 – 22 Março, 2002. Varricchio, S. L.; Gomes Jr., S.; Martins, N.; “s-Domain Approach to Reduce Harmonic Voltage Distortions Using Sensitivity Analysis”, IEEE Proceedings of the Winter Power Meeting, Columbus, Ohio, USA, 28 Janeiro – 1 Fevereiro, 2000. Gomes Jr., S.; Martins, N.; Portela C.; “Modal Analysis Applied to s-Domain Models of ac Networks”, IEEE PES Winter Meeting, Columbus, Ohio, Janeiro 2000. Varricchio, S. L.; Gomes Jr., S.; Martins, N.; “Two Powerful Network Modeling Approaches for the Modal Analysis of Harmonic Problems”, – VIII SEPOPE, Brasília, Brasil, 19 23 Maio, 2002. Lima, L. T. G.; Martins N.; Carneiro Jr., S.; “Dynamic Equivalents for Electromagnetic Transient Analysis Including Frequency-Dependent Transmission Line Parameters”, Proceedings of the IPST’97 International Power System Transients Conference, Seattle, USA, July, 1997. Varricchio, S. L.; Martins, N.; Lima, L. T. G.; Carneiro Jr., S.; “Studying Harmonic Problems Using a Descriptor System Approach”, Proceedings of the IPST’99 - International Conference on Power System Transients, Budapest, Hungary, June, 1999. Varricchio S. L.; Martins, N.; Lima, L. T. G.; “A NewtonRaphson Method Based on Eigenvalue Sensitivities to Improve Harmonic Voltage Performance”, IEEE Transactions on Power Delivery, vol. 18, no. 1, January 2003, pp. 334-342. Oualline, S.; “Practical C++ Programming”, O’Reilly & Associates, Inc., USA, 1995. Kruglinski J.D., (1997). “Inside Visual C++”, Redmond, Washington, Microsoft Press, 4a edition, 1997 Montenegro, F.; Pacheco, R.; “Orientação a Objetos em C++”, Editora Ciência Moderna, Brazil, 1994. Agostini, M. N.; Decker, I. C.; Silva, A. S.; “Desenvolvimento e Implementação de uma Base Computacional Orientada a Objetos para Aplicações em Sistemas de Energia Elétrica”, XIII Congresso Brasileiro de Automática, Florianópolis, Brasil, 2000. Manzoni, A.; Silva, A. S.; Decker, I. C.; “Power Systems Dynamic Simulation Using Object-Oriented Programming”, IEEE Transactions on Power Systems, Vol. 14, No. 1, February 1999, pp. 249−255. Araujo, L. R.; Garcia, P. A. N.; S., Pereira, J. L. R.; Carneiro Jr., S.; Vinagre, M. P.; “Modelagem Orientada a Objetos Aplicada na Solução de Sistemas de Distribuição”, XIV Congresso Brasileiro de Automática, Natal, Brasil, 2002. Chappel, D. “ActiveX and OLE”, Microsoft Press, USA, 1996. Rogerson, D.; “Inside COM”, Microsoft Press, USA, 1997. Task Force on Harmonics Modeling and Simulation, “Tutorial on Harmonics Modeling and Simulation”, IEEE Special Publication 98TP125-0, 1998. Alves, F. R. M.; Guimarães., C.; Pinto, H. J. P.; “Manual do Usuário do Programa ANAREDE V. 07”, CEPEL Relatório Técnico DPP/POL 530/87, Janeiro 1999. Diniz, R.; Gomes Jr., S.; Ferraz, J. C.; “Manual do Usuário do Programa ANATEM”, CEPEL, Janeiro 2002.
6
