Universidade Estadual de Santa Cruz - UESC
Departamento de Ciências Exatas e Tecnológicas - DCET
Programa de Pós-graduação em Física
MÉTODO ESPECTRONODAL DE GRADES
COMPOSTAS PARA PROBLEMAS DE AUTOVALOR
DE DIFUSÃO EM GEOMETRIA X,Y
‡
HUGO MENEZES DO NASCIMENTO VASCONCELOS
ILHÉUS, BA, BRASIL
Março de 2012
‡
Trabalho nanciado pela CAPES.
HUGO MENEZES DO NASCIMENTO VASCONCELOS
MÉTODO ESPECTRONODAL DE GRADES COMPOSTAS
PARA PROBLEMAS DE AUTOVALOR DE DIFUSÃO EM
GEOMETRIA X,Y
Dissertação apresentada ao Programa de PósGraduação em Física, Universidade Estadual
de Santa Cruz, para a obtenção do grau de
Mestre em Física.
Área de Concentração:
Ilhéus-BA, Brasil
2011
Física
Hugo Menezes do Nascimento Vasconcelos
A comissão examinadora, abaixo assinada, aprova a Defesa
MÉTODO ESPECTRONODAL DE GRADES COMPOSTAS
PARA PROBLEMAS DE AUTOVALOR DE DIFUSÃO EM
GEOMETRIA X,Y
Comissão Julgadora:
Dr. Dany S. Dominguez (UESC)
Dr. Francisco B. S. Oliveira (UESC)
Dr. Ricardo C. de Barros (UERJ)
Ilhéus-BA, 02 de março de 2012
Dedico este trabalho a minha família.
AGRADECIMENTOS
Agradeço a Deus por me permitir chegar até aqui;
Ao professor Dany Sanchez Dominguez pela paciência e pelos ensinamentos que vou
levar para toda vida, pois não foi fácil;
A mínha família pela força, que sempre me fez suportar todos os obstáculos;
Aos colegas do mestrado, em especial Raquel e Mirta, que sempre me zeram companhia no almoço e também nos pensamentos positivos para o futuro;
Aos amigos de longa data, que mesmo longe ou perto estavam sempre presentes: Alex,
Vivaldo, Leonardo, Joselito. E os novos amigos: Guilherme, Allen, Rafael, Daniel, que zeram
desse tempo mais interessante e não me deixaram desanimar;
Aos funcionarios da UESC, principalmente os do NBCGIB, PROFÍSICA, CPqCTR,
que sempre estavam disponíveis à ajudar;
A Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), que nanciou meus estudos.
Enm, agradeço a todos que direta ou indiretamente, contribuiram para minha chegada
até aqui.
NOMENCLATURA
J x,y
Corrente neutrônica na direção espacial
x
ou
Σa
Seção de choque macroscópica de absorção.
Σf
Seção de choque macroscópica de ssão.
ϕ
D
kef
ν
y.
Fluxo escalar de nêutrons.
Coeciente de difusão.
Fator de multiplicação efetivo.
Número médio de nêutrons liberado por evento de ssão
RESUMO
Nesta dissertação, o Método Espectronodal de Grades Múltiplas (CSG-SND-CN) é
desenvolvido com a nalidade de reduzir os custos algébricos e computacionais, dos métodos espectronodais existentes. Propomos o método espectronodal de grades compostas para
resolver numericamente problemas de autovalor de difusão em geometria
X, Y
(CSG-SND).
Este método utiliza uma grade espacial retangular grossa com uma célula por cada região do
problema. Primeiramente integramos transversalmente a equação de difusão, separadamente
em X e Y no interior de cada célula, e depois introduzimos aproximações constantes para
os termos de fuga transversal. Em seguida, usamos uma grade espacial na para discretizar
cada problema de difusão nodal unidimensional, que resolvemos separadamente.
Como as
direções espaciais são acopladas pelos termos de fuga transversal, utilizamos uma técnica de
direções alternadas para convergir a solução numérica. Assim, apresentamos nossos resultados
de modo a ilustrar a precisão e eciência do método aqui desenvolvido.
Palavras-chaves: reator nuclear, equação de difusão, problema de autovalor, espectronodal,
grades compostas
ABSTRACT
In this dissertation, a Composite Spatial Grid-Spectral Nodal Diusion - Constant
Nodal (CSG-SND-CN) is developed in order to reduce algebraic and computational costs
of existing spectral nodal methods.
We propose this method to numerically solve
X, Y
-
geometry one-speed diusion eigenvalue problems (CSG-SND). This method uses a coarse
rectangular spatial grid composed of one cell in each region of the problem. We rst integrate
transversely the diusion equation, separately in X and Y within each cell, and then we
introduce constant approximations for the transverse leakage terms. Next we use a ne spatial
grid to discretize each one-dimensional nodal diusion problem, which we solve separately.
As the spatial coordenate directions are coupled by the transverse leakage terms, we use a
technique of alternating directions to converge the numerical solution. We present numerical
results to illustrate the accuracy and eciency of the oened method.
Key-words: nuclear reactor, diusion equation, eigenvalue problem, spectral nodal, composite
grids
LISTA DE FIGURAS
2.1
Representação da lei de Fick.
2.2
Volume arbitrário onde migra uma população de nêutrons.
2.3
Condição de contorno e distância de extrapolação na superfície.
3.1
Domínio de cálculo que representa um núcleo retangular de dimensões
3.2
Grade espacial utilizada na discretização do domínio de cálculo.
3.3
Representação de um nodo
3.4
Representação das grandezas associadas a um nodo arbitrário
3.5
Grade espacial e incógnitas, considerando uma malha de
Γi,j
. . . . . . . . . . . . . . . . . . . . . . . . . . .
da grade espacial.
. . . . . . . . . . .
20
21
. . . . . . . .
25
Lx
Ly .
29
. . . . . . . .
30
. . . . . . . . . . . . . . . .
31
Γij
3 × 3 ; ( •)
e
. . . . . . .
33
grandezas
médias na célula espacial 3 por nodo; (◦) grandezas em arestas verticais 2 por
aresta; e (×) grandezas em arestas horizontais 2 por aresta. . . . . . . . . . . .
3.6
Aproximação do método DD para o uxo e a corrente no interior do nodo, a
esquerda discretização de malha grossa, a direita discretização de malha na. .
3.7
Representação das incógnitas do problema nas arestas verticais para uma linha
3.8
Representação das incógnitas do problema nas arestas horizontais para uma
linha
i.
Domínio da solução e grade espacial grossa exterior
4.2
Grades internas nas na direção
grossa
Λij .
x (Γx ),
Λ
e na direção
. . . . . . . . . . . . . .
y (Γy ),
Incógnitas dos sistemas (a) Problema unidimensional dependente de
y.
x,
. . . . . . . . . . . . . . . . . . . .
Relação entre as grades espaciais e o desvio relativo para um reator innito.
5.2
Desvio relativo percentual no cálculo do autovalor dominante para o reator
53
57
62
69
homogêneo nito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
72
Número de iterações externas segundo a grade espacial para o reator homogêneo
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
73
Geometria e condições de contorno para o problema de reator heterogêneo
simétrico.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
74
Desvio relativo percentual no cálculo do autovalor dominante para o problema
modelo do reator heterogêneo simétrico.
5.6
37
.
nito.
5.5
36
(b)
5.1
5.4
35
para uma célula
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Problema unidimensional dependente de
5.3
j.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1
4.3
34
. . . . . . . . . . . . . . . . . . . . .
76
Numero de iterações externas para o problema modelo do reator heterogêneo
simétrico.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
77
5.7
Geometria e condições de contorno para o problema do reator heterogêneo de
quatro regiões. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.8
Desvio relativo percentual no cálculo do autovalor dominante para o problema
modelo do reator heterogêneo de quatro regiões. . . . . . . . . . . . . . . . . .
5.9
78
80
Número de iterações externas para o problema modelo do reator heterogêneo
de quatro regiões. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
81
LISTA DE TABELAS
α
2.1
Coeciente Albedo
para as diferentes condições de contorno. . . . . . . . . .
3.1
Balanço de equações e incógnitas para os problemas unidimensionais de EDOs
gerados pela integração transversal. . . . . . . . . . . . . . . . . . . . . . . . .
5.1
41
Fator de multiplicação efetivo e número de iterações externas do problema
modelo de reator innito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2
26
68
Fator de multiplicação efetivo e número de iterações externas para o problema
modelo do reator homogêneo nito. . . . . . . . . . . . . . . . . . . . . . . . .
71
5.3
Propriedades físicas para o problema de reator heterogêneo simétrico.
74
5.4
Fator de multiplicação efetivo e numero de iterações externas para o problema
heterogêneo simétrico.
. . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.5
Propriedades físicas para o problema do reator de quatro regiões.
. . . . . . .
5.6
Fator de multiplicação efetivo e número de iterações externas para o problema
heterogêneo de quatro regiões. . . . . . . . . . . . . . . . . . . . . . . . . . . .
75
77
79
SUMÁRIO
1 INTRODUÇÃO
12
2 TEORIA DA DIFUSÃO
17
2.1
Taxa de Interação, Fluxo e Corrente de Nêutrons
. . . . . . . . . . . . . . . .
17
2.2
Lei de Fick
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
2.3
A Equação de Continuidade
2.4
Equação de Difusão e Condicões de Contorno
2.5
Limitações da teoria da difusão
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
. . . . . . . . . . . . . . . . . .
23
. . . . . . . . . . . . . . . . . . . . . . . . . .
26
3 MÉTODOS NUMÉRICOS PARA PROBLEMAS DE AUTOVALOR DE
DIFUSÃO
28
3.1
Equação de Difusão em Geometria X,Y . . . . . . . . . . . . . . . . . . . . . .
28
3.2
Equações de Balanço espacial
. . . . . . . . . . . . . . . . . . . . . . . . . . .
31
3.3
Método Diamond Dierence (DD) . . . . . . . . . . . . . . . . . . . . . . . . .
34
3.4
Método Espectro Nodal de Difusão com Aproximação Constante (SND-CN)
39
4 MÉTODO ESPECTRONODAL DE GRADES COMPOSTAS
52
4.1
Bases Matemáticas do Método CSG-SND-CN
. . . . . . . . . . . . . . . . . .
53
4.2
Esquema Numérico para o Método CSG-SND-CN . . . . . . . . . . . . . . . .
59
5 RESULTADOS NUMÉRICOS E DISCUSSÃO
67
5.1
Reator Homogêneo Innito . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
67
5.2
Reator Homogêneo Finito
69
5.3
Reator Heterogêneo Simétrico
. . . . . . . . . . . . . . . . . . . . . . . . . . .
73
5.4
Reator Heterogêneo de Quatro Regiões . . . . . . . . . . . . . . . . . . . . . .
77
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6 CONCLUSÕES E TRABALHOS FUTUROS
82
1 INTRODUÇÃO
A energia é um ingrediente essencial à vida humana desde as sociedades primitivas,
onde seu custo era praticamente zero por ser obtida de lenha das orestas. Porém, o consumo
energético foi crescendo de tal modo que outras formas de energia se tornaram necessárias.
Durante a Idade Média, rodas d'água começaram a ser usadas para tentar suprir as necessidades energéticas da população crescente. Depois da Revolução Industrial, a utilização de
carvão, petróleo e gás tornou-se necessária, mesmo tendo custos elevados para produção e
transporte [Goldemberg e Lucon 2007]. Nas últimas décadas a produção de energia através
de combustíveis fósseis consolidou-se como a principal alternativa energética, o que derivou
em graves problemas ambientais.
A geração de energia baseada em combustíveis fósseis libera na atmosfera grandes
quantidades de dióxido de carbono (CO2 ), este gás é o principal responsável pelo efeito estufa
[Mendonça e Gutierez 2000] e pelas mudanças climáticas provocadas pelo aquecimento global.
Nos dias atuais a humanidade enfrenta um grande desao: produzir energia para manter os
níveis de consumo, do estilo de vida moderno, e ao mesmo tempo diminuir o impacto das
atividades humanas no meio ambiente.
Produzir energia livre da emissão de
CO2
é uma
questão de sobrevivência.
Nos últimos anos, um grande número de especialistas e pesquisadores de renome da
área de energia nuclear a enxergam como uma fonte poderosa, que pode ser utilizada para
suprir as necessidades energéticas futuras, com baixo impacto ambiental.
Um dos motivos
principais para este posicionamento é que a energia nuclear não emite gases de efeito estufa
na atmosfera. Entretanto existem quatro grandes entraves que devem ser superados para a
expansão do uso da energia nuclear. Estes entraves são:
1. diminuir ou minimizar os estoques de materiais nucleares, que possam ser utilizados
com ns bélicos;
2. resolver o problema do armazenamento dos rejeitos radioativos de longa vida;
3. melhorar os indicadores econômicos da energia nuclear, principalmente os altos custos
de construção de novas usinas;
4. aprimorar a segurança e conabilidade operacional das usinas nucleares em operação ou
por construir, visando a possibilidade de excluir a ocorrência de acidentes com liberação
de material radioativo ao meio ambiente.
13
Temos que a produção de energia elétrica em reatores nucleares é sustentada principalmente pela reação de ssão nuclear de isótopos de urânio. Nesta reação um nêutron colide
com um núcleo pesado e provoca a ruptura deste núcleo em dois núcleos mais leves e simultaneamente são emitidos, em média, dois novos nêutrons. A energia liberada na reação, em
forma de energia cinética dos fragmentos de ssão, pode ser convertida em energia elétrica.
Para garantir a segurança e integridade do reator nuclear, a energia deve ser extraída do núcleo ao mesmo tempo em que as grandezas operacionais como temperatura e pressão devem
ser mantidas dentro dos níveis estabelecidos pelo projeto da instalação. Assim, os sistemas
de controle, de proteção e o comportamento dos componentes precisam ser cuidadosamente
testados e vericados, tanto em condições normais de operação, quanto em condições adversas
(prevendo acidentes ou o mau funcionamento).
Com esta nalidade, a física de reatores constitui uma área de pesquisa que entre
outras aplicações, desenvolve métodos de análise para predizer a estabilidade do reator, e o
seu comportamento durante diferentes tipos de transientes, inclusive em estado estacionário.
Uma parte destas análises são os cálculos de criticalidade, onde devemos computar o fator de
multiplicação efetivo à distribuição do uxo neutrônico. A distribuição de uxo é proporcional
a potência térmica liberada no núcleo do reator. Geralmente o problema associado a análises
de criticalidade é modelado matematicamente como um problema de autovalor, onde o fator
de multiplicação efetivo é denido como o autovalor dominante, e a distribuição de uxo é
a autofunção associada ao mesmo. Os cálculos de criticalidade são elementos essenciais para
garantir a operação segura de um reator nuclear.
Desta forma, a presente pesquisa visa contribuir na solução de um dos atuais entraves
da energética nuclear: a segurança operacional. Especicamente no desenvolvimento de um
método numérico que permita realizar cálculos de criticalidade com uso eciente dos recursos
computacionais.
Este método poderá ser utilizado na modelagem e simulação de estados
estacionários do núcleo do reator.
A migração de partículas neutras em um meio material é descrita pelo modelo matemático de transporte de nêutrons.
Este modelo é representado por uma equação integro
diferencial parcial linear, a qual considera que os nêutrons interagem com o meio sem alterar sua estrutura e não interagem entre si [Lewis e Miller 1993].
A equação de transporte
representa, portanto, um balanço entre produção e remoção de nêutrons, sendo que, em sua
generalidade, ela é uma equação dependente de sete variáveis: três espaciais, duas angulares,
uma de energia e uma temporal. O modelo de transporte apresenta uma alta complexidade
e sua resolução requer um conjunto de simplicações e aproximações.
14
Uma destas aproximações é a teoria da difusão, onde não se considera a dependência
angular do fenômeno de transporte. A equação de difusão apresenta uma redução do número
de variáveis independentes de sete para cinco: três espaciais, uma de energia e uma temporal.
Classicamente, a variável de energia, que é uma variável não-negativa contínua, é tratada
de forma discretizada em G intervalos contíguos, denominados grupos de energia (g = 1:G),
sendo convencional considerar g crescente para energia decrescente. Este método de discretização é conhecido como método multigrupo. A teoria multigrupo transforma a equação da
difusão de nêutrons em um sistema de G equações acopladas pelos termos de fonte de ssão e
espalhamento. Do ponto de vista computacional, a aproximação da difusão é menos robusta
que o modelo de transporte. Matematicamente, esta aproximação gera restrições na dependência angular da densidade de nêutrons. Fisicamente estas restrições podem ser sintetizadas
em duas premissas fundamentais [Duderstadt e Hamilton 1976]:
(a) os processos de espalhamento devem ser dominantes nas interações dos nêutrons com o
meio, isto é, o meio material tem que se apresentar fracamente absorvedor de nêutrons;
(b) a uência dos nêutrons deve ocorrer longe de fontes localizadas e de grandes descontinuidades materiais, pois nestas regiões do domínio espera-se maior anisotropia no processo
de transporte.
Contudo, a teoria da difusão tem sido amplamente usada em análises de reatores nucleares e, em geral, fornece resultados mais satisfatórios do que se poderia esperar teoricamente.
Em algumas situações, em reatores nucleares de potência, onde ocorrem grandes gradientes, a
aproximação da difusão pode falhar, como por exemplo, em regiões dos reatores térmicos onde
as barras de controle estão inseridas, ou nas proximidades do contorno em reatores compactos.
Ademais, em razão da complexidade ou até da impossibilidade do tratamento analítico completo da equação da difusão utilizada para modelar a distribuição de nêutrons em
reatores nucleares, métodos numéricos são desenvolvidos e aplicados à equação da difusão em
determinada geometria. Esses métodos numéricos discretizam as variáveis espaciais, e usam
vários esquemas diretos e iterativos para resolver o sistema de equações algébricas resultante
conforme mencionamos acima, a variável energia é tratada pela aproximação multigrupo. A
variável espacial pode ser discretizada por métodos de malha na, como os métodos de diferenças nitas (DF) [Lamarsh e Baratta 2001]; métodos de malha média como os métodos de
elementos nitos [Zienkiewicz 1971] ou os métodos de malha grossa, entre os quais estão os
métodos nodais [Lawrence 1986] e [Badruzzaman 1990].
15
Por razões vinculadas à eciência computacional, ao discretizarmos as variáveis espaciais, é desejável utilizarmos métodos de malha grossa; isto é, utilizar grades espaciais com
células largas, o que conduz a um menor número de incógnitas do problema. Dentre os métodos de malha grossa destacam-se os Métodos Nodais (MN). Os MN são baseados no conceito
de integração transversal, onde as equações de difusão são integradas separadamente em cada
uma das direções espaciais. Em seguida, os termos de fonte e os termos de fuga transversal
são aproximados por polinômios de baixa ordem. Entre os métodos nodais, aparecem o método Nodal Constante (CN) [Lawrence 1986] onde as aproximações são polinômios de ordem
zero, e o método nodal linear (LN) onde as aproximações são polinômios de primeira ordem
[Azmy 1988].
Dentro da família dos métodos nodais, destacam-se os métodos espectronodais (SGF).
Nestes métodos, apenas os termos de fuga transversal são aproximados por polinômios de
baixa ordem, os termos de fonte são tratados analiticamente dentro da formulação. Os métodos espectronodais foram desenvolvidos a partir dos anos 90, iniciando com o método espectronodal constante (SGF-CN) na formulação de ordenadas discretas (Sn ) da equação de
transporte apresentado por [Barros e Larsen 1992] para problemas de fonte xa. Posteriormente [Barros e Silva 1999] apresentaram o método espectronodal constante (SGF-SD-CN)
para problemas
Sn
de autovalor. Mais recentemente os autores [Barros et al. 2003] apresen-
taram o método espectronodal constante (SND-CN) para problemas de autovalor de difusão.
O método espectro nodal linear (SGF-LN) para problemas
Sn
de fonte xa foi apresentado
por [Dominguez e Barros 2007].
De forma geral, os métodos espectronodais mostraram-se menos sensíveis aos erros de
discretização espacial que os métodos nodais convencionais; isto é, eles podem gerar resultados numéricos precisos para malhas mais grossas. No entanto, os métodos espectronodais
apresentaram duas grandes desvantagens: o alto custo algébrico para obter os esquemas iterativos de solução numérica, e a necessidade de implementar algoritmos computacionais mais
complexos. Visando diminuir o custo algébrico e computacional dos métodos espectronodais,
aparecem os métodos de grades compostas (CSG). Estes métodos aplicam o conceito de integração transversal em uma grade grossa, que deve coincidir com as regiões do domínio e
aproximam as fugas transversais nesta grade externa. Em seguida, os problemas "unidimensionais" resultantes são discretizados utilizando uma grade mais na, e resolvidos utilizando
formulações SGF convencionais. A discretização das variáveis espaciais é feita em duas etapas:
I. em nível de região material;
16
II. e no interior das regiões materiais.
Os métodos de grades compostas apresentaram menor custo algébrico e computacional
que os métodos espectronodais convencionais. Entretanto, enfatizamos que nos método CSG
as fugas são aproximadas na região, e esta aproximação é mais grosseira que a aproximação
dos métodos SGF. Mesmo assim, em [Dominguez et al. 2009] mostrou-se que os métodos
espectronodais de grades compostas são apropriados para calcular grandezas integrais em
problemas de fonte xa, formulação
Sn .
Consequentemente deve-se obter bons resultados no
cálculo do fator de multiplicação efetivo e dos valores médios de distribuição de potência em
problemas de autovalor de difusão.
Neste trabalho, propomos desenvolver o método espectronodal de grades compostas
para problemas de autovalor, utilizando a aproximação monoenergética da difusão em geometria cartesiana bidimensional. Para aproximar as fugas transversais nas regiões do domínio,
utilizaremos polinômios de ordem zero.
Esta nova formulação permitirá desenvolver roti-
nas computacionais, que deve obter resultados precisos no cálculo de grandezas integrais
em malhas grossas, oferecendo resultados com precisão similares aos métodos espectronodais
clássicos, com menor custo computacional.
A presente dissertação está dividida em cinco capítulos. No presente capítulo introdutório é oferecida uma justicativa da necessidade desta pesquisa e são descritos os objetivos
do trabalho.
No capítulo 2 mostramos os fundamentos físicos e o modelo matemático da
difusão de nêutrons.
No capítulo 3 são apresentados alguns métodos numéricos utilizados
para resolver problemas de autovalor de difusão; especicamente mostramos o método de
diferenças nitas, e o método espectronodal de difusão com aproximação constante. No capítulo 4 oferecemos os fundamentos matemáticos do método espectronodal de difusão de grades
múltiplas, objetivo principal deste trabalho. No capítulo 5 mostramos os resultados para problemas modelos típicos em domínios homogêneos, e para um reator heterogêneo simétrico e
de quatro regiões. Finalmente, no capítulo 6 comentamos os resultados obtidos neste trabalho
e oferecemos perspectivas para novos trabalhos.
2 TEORIA DA DIFUSÃO
Um dos principais problemas na física de reatores nucleares é a determinação da distribuição de nêutrons no interior do reator [Lamarsh e Baratta 2001]. A distribuição de nêutrons
determina a taxa de reação com que as diversas interações ocorrem. Ademais, o estudo do
comportamento da população de nêutrons permite caracterizar a estabilidade da reação de
ssão em cadeia. Para estimar o perl da distribuição de nêutrons dentro do reator devemos
pesquisar o processo de transporte de nêutrons, ou seja, o movimento dos nêutrons e seu uxo
no interior do reator.
Um nêutron que migra no núcleo do reator poderá sofrer interações
de espalhamento com os núcleos dos átomos do meio, que leve a ssão do núcleo ou não, ou
ainda escapar do núcleo do reator.
A maioria dos estudos sobre reatores nucleares tratam o movimento dos nêutrons
como um processo de difusão, assumindo que os nêutrons se difundem de regiões com alta
densidade de nêutrons para regiões com baixa densidade. Uma analogia utilizada pela maioria
dos autores é considerar o movimento dos nêutrons semelhante à difusão de um gás molecular
rarefeito, onde as moléculas do gás em movimento se difundem para minimizar as variações
na concentração. A teoria da difusão é limitada a condições especícas sobre o meio material
e a características do uxo. Entretanto, do ponto de vista matemático e conceitual é muito
mais simples que os modelos de transporte.
A teoria da difusão tem oferecido resultados
satisfatórios para problemas de física de reatores durante décadas, sendo o principal modelo
utilizado para o desenho e operação de reatores na atualidade.
No presente capítulo, apresentamos os fundamentos físicos e matemáticos da equação
de difusão de nêutrons.
Para começar, introduzimos as grandezas físicas de uxo, taxa de
reação e corrente de nêutrons. Em seguida apresentamos a lei de Fick e a equação de continuidade, que se combinam para conformar a equação de difusão. Por último, apresentamos as
condições de contorno da equação de difusão e alguns comentários sobre as limitações desta
teoria. Um leitor experiente em física dos reatores pode desconsiderar a leitura deste capítulo.
2.1 Taxa de Interação, Fluxo e Corrente de Nêutrons
Nosso objetivo é determinar a distribuição de nêutrons no núcleo do reator; para isso
devemos considerar o movimento dos nêutrons durante sua migração no núcleo e as interações
dos nêutrons com os núcleos atômicos do meio. Denimos a densidade de nêutrons
em um ponto
⃗r
n(⃗r, E, t)
do núcleo do reator como o número esperado de nêutrons no volume
dV
18
relativo a
⃗r
com energias entre
E
e
(E + dE)
no instante de tempo
t.
A densidade de nêutrons nos permite calcular a taxa de reação em qualquer ponto
do núcleo do reator. A frequência de interação é proporcional à probabilidade de ocorrência
de uma reação no meio material (seção de choque macroscópica total Σt )
e a velocidade
(energia) do nêutron
frequência de interação
denimos a taxa de reação
⃗r
no intervalo de energia
F
= Σt v
como o valor esperado da taxa de interação que ocorre na posição
[E, E + dE]
no instante de tempo
t
como
F (⃗r, E, t) = ΣT (⃗r, E)vn(⃗r, E, t)
o produto da velocidade do nêutron pela densidade de nêutrons em um ponto
(2.1)
⃗r
aparece
frequentemente em física dos reatores e é denido como uxo escalar de nêutrons
ϕ(⃗r, E, t) = vn(⃗r, E, t) ,
que é expresso nas unidades
[cm−2 s−1 ].
O uxo escalar de nêutrons pode ser interpretado
como o número total de nêutrons com energia
qualquer direção no instante
(2.2)
E
que atravessam uma unidade de área em
t.
Neste trabalho consideramos nêutrons monoenergéticos; isto é, todos os nêutrons
propagam-se à mesma velocidade; esta aproximação permite desconsiderar a dependência
energética do uxo e das taxas de reação. Obtemos o uxo de nêutrons em aproximação de
uma velocidade integrando a eq. (2.2) em todo o domínio possível de energia
∫
∞
ϕ(⃗r, t) =
ϕ(⃗r, E, t)dE,
(2.3)
0
nalmente podemos escrever a taxa de reação total para nêutrons monoenergéticos na forma
F (⃗r, t) = ΣT (⃗r)ϕ(⃗r, t).
Aqui, a eq.
(2.4)
(2.4) refere-se à taxa de interação total; isto é, inclui todos os tipos
de reações nucleares nêutronnúcleo. Taxa de interações para reações especicas, possuem
expressões semelhantes. Assim, o número de colisões de espalhamento no ponto
calculado como
⃗r
pode ser
19
Fs (⃗r, t) = Σs (⃗r)ϕ(⃗r, t) ;
e o número de nêutrons absorvidos no volume
dV
no ponto
(2.5)
⃗r
pode ser expresso na forma
Fa (⃗r, t) = Σa (⃗r)ϕ(⃗r, t) .
(2.6)
Maiores detalhes sobre as grandezas envolvidas na teoria da difusão podem ser encontrados em [Duderstadt e Hamilton 1976] e [Lamarsh e Baratta 2001].
2.2 Lei de Fick
Estando o uxo e a corrente neutrônica estreitamente relacionados, devido à complexidade do fenômeno de transporte de nêutrons, não é uma tarefa simples descrever matematicamente a relação entre estas duas grandezas.
A teoria da difusão estabelece uma relação entre uxo e corrente através da lei de Fick,
que foi originalmente utilizada para descrever o processo de difusão química. Em uma solução
química se a concentração de soluto é maior em uma região que em outra, o soluto difunde-se
da região de maior concentração para a região de menor concentração. Além disso, vericouse que a taxa de uxo de soluto é proporcional ao negativo do gradiente da concentração de
solutos. Esta é a declaração original da lei de Fick. A principal aproximação da teoria da
difusão estabelece que os nêutrons que migram no interior do núcleo do reator comportam-se
da mesma forma que o soluto em uma solução química. Assim, se a densidade de nêutrons é
maior em uma parte de um reator que em outra, há uma corrente de nêutrons na direção da
região de menor densidade de nêutrons. Por exemplo, suponha que o uxo varia ao longo da
direção espacial
x,
como mostrado na gura (2.1) [Lamarsh e Baratta 2001]. Então a lei de
Fick na teoria da difusão pode ser escrita na forma
Jx (x, t) = −D(x)
Nesta expressão,
Jx
dϕ(x, t)
.
dx
(2.7)
é igual ao número líquido de nêutrons que passam por unidade de
tempo através de uma unidade de área perpendicular à direção
x.
O parâmetro
equação (2.7) é denominado de coeciente de difusão e tem unidades de
cm.
D
na
O coeci-
ente de difusão é uma aproximação para ajustar a teoria da difusão à teoria de transporte
[Duderstadt e Hamilton 1976].
A eq. (2.7) mostra que, se existe um gradiente de uxo negativo na direção positiva
20
Figura 2.1: Representação da lei de Fick.
do eixo
x,
então o uxo diminui nessa direção espacial. Em outras palavras, há uma corrente
de nêutrons ao longo da direção
x
positiva, tal como indicado na gura 2.1. Para entender a
origem desta corrente, considere os nêutrons que passam através do plano perpendicular ao
eixo das abscissas em
x = 0.
Existem nêutrons que passam através do plano da esquerda para
a direita (no sentido positivo do eixo
(no sentido negativo do eixo
x).
x) e nêutrons que atravessam da direita para a esquerda
No entanto, uma vez que a densidade de nêutrons e o uxo
são maiores para valores negativos de
do que à direita do plano
x = 0.
x,
há mais colisões por unidade de volume à esquerda
Portanto, mais nêutrons são espalhados da esquerda para
a direita do que o inverso. Isto acarreta que há uma corrente líquida de nêutrons na direção
positiva do eixo
x,
tal como previsto pela eq. (2.7).
No caso mais geral, a posição
⃗r
é descrita por três variáveis espaciais, e consequente-
mente, o uxo é geralmente uma função de três variáveis espaciais. No caso geral, a lei de
Fick aparece na forma
⃗ r, t) = −D(⃗r)∇Φ(⃗r, t)
J(⃗
onde
∇
(2.8)
é o operador gradiente na geometria do problema.
2.3 A Equação de Continuidade
Considere
V
um volume arbitrário conforme a gura (2.2), este volume é formado pelos
pontos interiores do volume
Vi (⃗r)
e pela superfície
Γ(⃗r)
que limita o contorno do volume. No
interior do volume migra uma população de nêutrons e a variação no tempo da densidade
de nêutrons no interior de
conforme a expressão
Vi
é inuenciada por fenômenos de perda e produção de nêutrons
21
Figura 2.2: Volume arbitrário onde migra uma população de nêutrons.
[Taxa de variação do número de nêutrons em V]
= [Produção]
- [Perda] .
Os processos associados à remoção de nêutrons do interior do volume, são reações nucleares de
absorção e a fuga de nêutrons do domínio, que acontece quando nêutrons atravessam para fora
da superfície
Γ(⃗r).
Os processos associados à produção de nêutrons são a existência de fontes
xas no interior do domínio e as fontes de ssão, se o meio for multiplicativo. Inicialmente
agrupamos estes fenômenos de produção de nêutrons, com o nome de fonte, sem caracterizar
a origem dos nêutrons. Considerando os fenômenos descritos podemos considerar a seguinte
equação de balanço
[Taxa de variação do número de nêutrons em V]=
[Taxa de produção de nêutrons em V]
−
[Taxa de absorção de nêutrons em V]
−
fugas de nêutrons em V] .
[Taxa de
(2.9)
A seguir, formulamos matematicamente cada um dos termos que aparecem na eq. (2.9). Seja
n(⃗r, t)
a densidade de nêutrons em qualquer ponto do volume
número total de nêutrons em
V
V
no instante de tempo
t.
O
é calculado como
∫
n(⃗r, t)dV
V
onde a integral de volume indica que a integração deve ser realizada por todo
variação do número de nêutrons no interior de
V
é
V.
A taxa de
22
∫
Taxa de variação do número de nêutrons em V
=
V
Em seguida, denimos
cm3
V
a partir de fontes em
V.
s(⃗r, t)
∂n(⃗r, t)
dV .
∂t
como sendo a taxa na qual os nêutrons são emitidos por
A taxa com que os nêutrons são produzidos em todo o volume
é dada por
∫
=
Taxa de Produção
s(⃗r, t)dV ;
V
Em continuidade, a taxa na qual os nêutrons são perdidos por reações de absorção por
cm3 /s
é igual à taxa de reações de absorção denida na eq. (2.6). Se considerarmos todo o volume
V,
a perda total de nêutrons por segundo, devido à absorção é dada por:
∫
Taxa de Absorção
∫
=
Σa (⃗r)ϕ(⃗r, t)dV =
V
Σa (⃗r)vn(⃗r, t)dV .
V
Para avaliarmos as fugas, considere o uxo de nêutrons nas proximidades da superfície de
contorno
Γ
que circunda o volume
na superfície
Γ
e
⃗n
V.
Seja
J(⃗r, t)
o vetor densidade de corrente de nêutrons
um vetor unitário normal apontando para fora da superfície e, de acordo
com os resultados da seção anterior podemos armar que
⃗ r, t) · ⃗n
J(⃗
é o número líquido de nêutrons que atravessa a superfície do interior do domínio ao exterior
por unidade de área
(cm2 /s).
Daqui resulta que, a taxa total de fuga de nêutrons (que pode
ser positiva ou negativa) através da superfície
Γ
do volume é
∫
Taxa de Fuga
⃗ r, t) · ⃗ndΓ ,
J(⃗
=
Γ
esta integral de superfície pode ser transformada em uma integral de volume usando o teorema
da divergência de Gauss. Assim,
∫
∫
⃗ r, t) · ⃗ndΓ =
J(⃗
⃗ r, t)dV .
div J(⃗
Γ
V
Então,
∫
Taxa de Fuga
=
divJ(⃗r, t)dV .
V
Finalmente, podemos obter a equação da continuidade através da introdução dos resultados
23
anteriores para cada um dos termos da eq. (2.9). Após a substituição, o resultado aparece na
forma
∫
V
∂n(⃗r, t)
=
∂t
∫
∫
s(⃗r, t)dV −
V
∫
Σa (⃗r)ϕ(⃗r, t)dV −
V
⃗ r, t)dV .
div J(⃗
V
Todas as integrais acima devem ser executadas sobre o mesmo volume, ademais, a
equação deve ser válida em qualquer volume arbitrário. Assim, a igualdade deve ser satisfeita
também para os integrandos. O resultado nal aparece como
∂n(⃗r, t)
⃗ r, t) .
= s(⃗r, t) − Σa (⃗r)ϕ(⃗r, t) − div J(⃗
∂t
(2.10)
A eq. (2.10) é a forma geral da equação monoenergética da continuidade dependente do tempo,
representando o balanço de nêutrons no interior de um volume arbitrário. Para problemas
estacionários ou de equilíbrio, onde a densidade de nêutrons não é uma função do tempo, esta
equação se reduz à
⃗ r) + Σa (⃗r)ϕ(⃗r) = s(⃗r) ,
div J(⃗
(2.11)
que é conhecida como a equação monoenergética de estado estacionário de continuidade. Um
elemento importante na equação anterior é o termo de fonte
s(⃗r, t)
que no caso mais geral
pode ser escrito na forma
s(⃗r) =
1
νΣf (⃗r)ϕ(⃗r) + Q(⃗r) ,
kef
(2.12)
onde o primeiro termo representa a fonte de ssão e o segundo a fonte externa de nêutrons.
2.4 Equação de Difusão e Condicões de Contorno
A continuidade para problemas estacionários é dada pela eq. (2.11), a mesma caracteriza o balanço neutrônico em um volume arbitrário e pode ser utilizada como modelagem
matemática do problema de transporte de nêutrons monoenergéticos. Infelizmente, a equação
da continuidade tem duas incógnitas, o uxo de nêutrons,
ϕ(⃗r), e a corrente de nêutrons, J(⃗r).
Para obtermos solução única, precisamos de uma relação entre essas duas grandezas. Essa
relação foi mostrada na seção 2.2, sendo chamada de lei de Fick, eq. (2.8). A equação da
difusão para problemas estacionários pode ser expressa como um sistema de equações em derivadas parciais de primeira ordem formado pela equação de continuidade (2.11) e pela lei de
Fick (2.8); algumas formulações numéricas partem diretamente deste sistema de equações em
24
sua construção. Outra alternativa é substituirmos a lei de Fick na equação de continuidade,
para obter uma equação em derivadas parciais de segunda ordem. Esta equação é chamada de
equação da difusão de nêutrons, e supondo que o coeciente de difusão
D
não é uma função
da posição, o resultado aparece na forma
D∇2 ϕ(⃗r) − Σa (⃗r)ϕ(⃗r) =
onde o símbolo
∇2 = div · grad
1
νΣf (⃗r)ϕ(⃗r) + Q(⃗r) ,
kef
(2.13)
é chamado de laplaciano, que envolve derivadas de segunda
ordem e depende da geometria do problema. Em volumes onde não existem fontes externas
de nêutrons, como exemplo, no núcleo de reatores de potência, o último termo da eq. (2.13)
é nulo.
O uxo de nêutrons pode ser encontrado resolvendo a equação de difusão.
Como a
equação de difusão é uma equação diferencial parcial, para obtermos solução única é necessário especicar certas condições adicionais que devem ser satisfeitas pela solução. Algumas
destas condições são determinadas a partir de requisitos óbvios, para obtermos um uxo físico
razoável.
Por exemplo, desde que um uxo negativo ou imaginário não tem sentido físico,
conclui-se que
ϕ
deve ser uma função real e não negativa. O uxo também deve ser nito,
exceto talvez nos pontos singulares articiais de uma distribuição de fonte.
Em muitos problemas de migração de nêutrons, existe uma superfície externa que
separa o meio material e o ambiente. A lei de Fick não é válida nas imediações desta superfície,
e consequetemente a equação de difusão não é válida. No entanto, métodos de ordem superior
mostram que o uxo calculado a partir da equação de difusão desaparece a uma pequena
distância
d
além da superfície. Então, o uxo determinado a partir da equação de difusão é
quase igual ao uxo exato no interior do meio [Lamarsh e Baratta 2001]. A suposição de que
o uxo desaparece a uma pequena distância
d
para além da superfície não tem sentido físico.
Na verdade, é uma aproximação matemática conveniente, que fornece um grau de precisão
aceitável das estimativas do uxo. Esta situação é ilustrada na gura (2.3). O parâmetro
d
é
conhecido como a distância de extrapolação, e na maioria dos casos de interesse é dada pela
fórmula simples
d = 0, 71λtr ,
onde
λtr
é o livre caminho médio de transporte do meio; isto é
λtr = 3D ,
(2.14)
25
Figura 2.3: Condição de contorno e distância de extrapolação na superfície.
e assim
d
torna-se:
d = 2, 31D .
A maioria dos valores medidos do coeciente de difusão
(2.15)
D, para materiais não gasosos,
são geralmente menores que um centímetro. Assim, da eq. (2.15), vê-se que a distância de
extrapolação
d,
é geralmente muito pequena em comparação com as dimensões do reator.
Portanto, muitas vezes, é possível ao resolver a equação de difusão, assumir que o uxo
desaparece na superfície real do sistema.
Nos casos em que
d
não é desprezível, é preciso
estabelecer um limite matemático para o problema; esse limite é conhecido como o limite
extrapolado, e está localizado a uma distância
d
a partir do limite físico real. Essa condição
de contorno é conhecida como a condição de contorno de vácuo. É normalmente aplicada em
problemas onde existe um meio material com densidade baixa [Lamarsh e Baratta 2001].
Uma alternativa para expressar as condições de contorno é considerarmos condições
de contorno do tipo Albedo.
Nestas condições de contorno faz-se, a corrente neutrônica
(proporcional ao gradiente de uxo) na fronteira ser proporcional ao uxo de nêutrons no
contorno.
O coeciente de proporcionalidade
α
é conhecido como coeciente Albedo; as
condições de contorno Albedo são expressas na forma
J(Γ) = αϕ(Γ) ,
(2.16)
os coecientes Albedo para as diferentes condições de contorno são mostrados na tabela 2.1.
A condição de contorno reexiva aparece em fronteiras onde não existe gradiente de uxo; isto
acontece em problemas com condições de simetria. A condição de contorno de vácuo aparece
quando a densidade do meio, do outro lado do contorno, é muito menor que a densidade do
26
núcleo do reator, esta condição suporta a anulação do uxo na distância extrapolada, conforme
discutido acima. A condição de contorno de uxo zero está associada a núcleos de reatores
muito grandes, onde a distância de extrapolação pode ser desconsiderada; e a condição de
contorno de reexão abundante deve ser utilizada em reatores com reetor, onde o tamanho
do reetor deverá ser maior que vários livres caminhos médios.
Maiores detalhes sobre as
condições de contorno podem ser encontradas em [Stacey 2001] e [Barros e Filho 2007].
Tabela 2.1: Coeciente Albedo
α
para as diferentes condições de contorno.
Tipo
α
Condição de Contorno reexiva
0
Condição de Contorno tipo vácuo
0.5
Condição do tipo uxo zero no contorno
Condição de reetor W abundante
√ ∞
DW ΣaW
Também é necessário especicar condições adicionais em uma interface entre dois diferentes meios de difusão, como a interface entre o elemento combustível e o canal de moderação.
Uma vez que os nêutrons atravessam a interface sem sofrer colisões, não é difícil perceber que
tanto o uxo quanto a componente da corrente neutrônica normal à superfície, devem ser contínuos em toda a interface. Assim, em uma interface entre duas regiões
A
e
B,
as seguintes
relações devem ser satisfeitas:
onde
ϕA
(JB )n
e
ϕB
ϕA = ϕB
(2.17)
(JA )n = (JB )n
(2.18)
são, respectivamente, o uxo nas regiões
AeB
avaliados na interface, e
(JA )n
e
são componentes normal da corrente de nêutrons avaliadas na interface. As eq. (2.17)
e (2.18) são denidas como condições de continuidade na interface.
2.5 Limitações da teoria da difusão
A Lei de Fick expressa o fato de que, se o gradiente do uxo é negativo, então a
densidade da corrente é positiva.
Isto signica que as partículas se difundem da região de
maior uxo para a região de menor uxo, tendo colisões no meio.
A Lei de Fick impõe
algumas limitações sob as seguintes condições:
• Proximidade com as fronteiras:
A derivação assumiu um meio innito.
Para um meio nito, a lei de Fick é válida
27
somente em pontos que são maiores do que alguns livres caminhos médios, a partir das
bordas do meio. Isto é assim, uma vez que o termo exponencial decresce rapidamente
com a distância, tendendo a zero, e os pontos de apenas alguns livres caminhos médios,
a partir do ponto onde o uxo é calculado, fazem uma contribuição signicativa para a
integral.
• Proximidade com as fontes ou dissipadores
Foi assumido que a contribuição para o uxo é dada principalmente, a partir de colisões
de espalhamento. Fontes podem estar presentes, por causa do fator de atenuação, no
entanto, um pequeno número de nêutrons de fontes de origem contribuirá para o uxo,
se forem maiores do que alguns livres caminhos médios.
• Sistema LAB de espalhamento anisotrópico
O espalhamento isotrópico no sistema LAB deve ocorrer a baixas energias, o que normalmente não acontece. Entretanto, a lei de Fick ainda é válida com anisotropia moderada
no espalhamento, se uma forma modicada do coeciente de difusão é utilizada, baseado
na teoria de transporte.
• Meios fortemente absorventes
O uxo deve variar suavemente, o que não ocorre em meios fortemente absorvedores.
Assim, a lei de Fick aplica-se aos sistemas em
Σa << Σs .
• Proximidade com a interface
A suposição de um meio uniforme foi usada na derivação da lei de Fick. Na fronteira
entre dois meios, de propriedades de espalhamento diferentes, a lei de Fick ainda é
válida, desde que a mudança brusca não conduza a um uxo que varia rapidamente;
o que invalida as expansões da série de Taylor do uxo na derivação. As derivadas de
terceira ordem, de fato contribuem para o integrante do
J.
Assim a lei de Fick seria
válida se a segunda derivada do uxo não variasse sensivelmente.
• Fluxo variando no tempo rapidamente
Foi assumido que o uxo é independente do tempo. É possível relaxar esta exigência,
se a mudança em
caminhos médios.
ϕ
é pequena durante o tempo que um nêutron viaja alguns livres
3 MÉTODOS NUMÉRICOS PARA PROBLEMAS DE
AUTOVALOR DE DIFUSÃO
Problemas de autovalor de difusão estão ligados a cálculos de criticalidade em reatores nucleares.
Estes cálculos são fundamentais para garantirem a estabilidade e segurança
operacional do reator nuclear. O modelo matemático geralmente utilizado é a equação estacionária da difusão que como foi mostrado no capítulo anterior é representada por um sistema
de equações diferenciais parciais (EDPs). A solução analítica da equação de difusão apenas
pode ser obtida em geometrias muito simples, geralmente carentes de valor para aplicações
práticas. Para problemas heterogêneos em geometrias que descrevam com precisão o núcleo
do reator, métodos numéricos devem ser utilizados.
Os métodos numéricos transformam o sistema de EDPs em um sistema de equações
lineares e algébricas (SELA) que pode ser resolvido utilizando métodos diretos ou iterativos
segundo as características do sistema. Em geral os métodos numéricos que resolvem problemas
de autovalor de difusão seguem um processo que compreende as seguintes etapas: (i) discretização do domínio de interesse; (ii) obtenção das equações discretizadas de balanço espacial;
(iii) introdução de aproximações para os operadores diferenciais, para os termos de fonte e de
fuga transversal, e equações auxiliares; (iv)
montagem do esquema numérico de varredura
e construção do SELA associado; e (v) resolução do problema de autovalor combinando o
método da potência com um método de resolução de SELAs.
No presente capítulo, mostramos as etapas seguidas na construção de um método
numérico para resolver problemas de autovalor de difusão em geometria cartesiana bidimensional. Especicamente, mostramos o método
Diamond Dierence (DD) e o método espectro
nodal de difusão com aproximação constante (SND-CN,
Spectral Nodal Diusion - Constant
Nodal).
3.1 Equação de Difusão em Geometria X,Y
Consideramos que nosso domínio de interesse compreende um meio multiplicativo retangular de comprimento
Lx
e altura
Ly
conforme ilustrado na gura 3.1, no interior deste
meio não existem fontes externas. A equação da difusão de nêutrons em geometria
aproximação de uma velocidade aparece na forma
X, Y
e
29
∂J x (x, y) ∂J y (x, y)
1
+
+ Σa (x, y)ϕ(x, y) =
νΣf (x, y)ϕ(x, y) ,
∂x
∂y
kef
∂ϕ(x, y)
, 0 ≤ x ≤ Lx
∂x
(3.2)
∂ϕ(x, y)
, 0 ≤ y ≤ Ly ,
∂y
(3.3)
J x (x, y) = −D(x, y)
J y (x, y) = −D(x, y)
onde
x
e
y
são as variáveis independentes do sistema representando a posição através das
direções espaciais,
espacial
ssão,
Ju (x, y)
u, Σa (x, y)
nêutrons,
(3.1)
kef
onde
u = (x
ou
y)
a seção de choque macroscópica de absorção,
o fator de multiplicação efetivo,
Σf (x, y)
representa a corrente neutrônica na direção
ν
ϕ(x, y)
uxo escalar de
o número provável de nêutrons por evento de
a seção de choque macroscópica de ssão e
D(x, y)
o coeciente de difusão.
Figura 3.1: Domínio de cálculo que representa um núcleo retangular de dimensões
Lx
e
Ly .
No sistema de equações de difusão consideramos condições de contorno do tipo albedo.
Temos uma condição de contorno para cada uma das fronteiras do domínio (vide gura 3.1); ou
seja, fronteira esquerda (L
Left), direita (RRight), inferior (B Bottom) e superior (T Top).
As condições de contorno são descritas pelas expressões
J˜x (0, y) = −αL ϕ̃(0, y) ,
(3.4)
J˜x (Lx , y) = αR ϕ̃(Lx , y) ,
(3.5)
Jˆy (x, 0) = −αB ϕ̂(x, 0) ,
(3.6)
30
Jˆy (x, Ly ) = αT ϕ̂(x, Ly ) .
(3.7)
O primeiro passo na construção de um método numérico é a discretização espacial do
domínio de cálculo. Para isso introduzimos uma grade espacial ou malha com
direção
x,
e
J
é formada por
comprimento
divisões na direção
I×J
hx i
nodos
e altura
Γi,j ,
hyj .
y,
I
divisões na
conforme ilustrado na gura 3.2. Esta grade espacial
com
i = 1 : I
e
j = 1 : J,
cada um destes nodos tem
As especicações geométricas do nodo
Γi,j
são mostradas na
gura 3.3.
Figura 3.2: Grade espacial utilizada na discretização do domínio de cálculo.
Como o núcleo do reator é composto por regiões com diferentes composições, as constantes físico-neutrônicas
(D, Σa e νΣf ) são descritas por funções contínuas por partes, a discre-
tização do problema deve respeitar a restrição que no interior de um nodo
Γij
as propriedades
do meio sejam constantes. Escrevemos então o sistema de equações de difusão, eq.(3.1, 3.2 e
3.3), para um nodo arbitrário
da grade espacial na forma
Γij = {(x, y) ∈ R2 | xi−1/2 ≤ x ≤ xi+1/2 e yj−1/2 ≤ y ≤ yj+1/2 }
31
Figura 3.3: Representação de um nodo
Γi,j
da grade espacial.
1
∂J x (x, y) ∂J y (x, y)
+
+ Σaij ϕ(x, y) =
νΣfij ϕ(x, y) ,
∂x
∂y
kef
∂ϕ
(x, y),
∂x
xi−1/2 ≤ x ≤ xi+1/2
∂ϕ
(x, y),
∂y
yj−1/2 ≤ y ≤ yj+1/2 .
J x (x, y) = −Dij
J y (x, y) = −Dij
As eq.
(3.8, 3.9 e 3.10) representam o sistema de equações de difusão para o nodo
que nos leva a
I ×J
(3.8)
(3.9)
(3.10)
Γij ,
o
sistemas de EDPs acoplados pelas condições de continuidade de uxo
e corrente nas interfaces dos nodos. Destacamos que mesmo com o aumento do número de
equações, com relação ao sistema descrito pelas Eqs. (3.1, 3.2 e 3.3), estes sistemas apresentam
coecientes constantes o que simplica a resolução dos mesmos.
Finalizamos a etapa de discretização do domínio com a discretização das condições de
contorno. Considerando uma grade espacial de
I ×J
nodos conforme apresentada nesta seção
as condições de contorno, eq. (3.4-3.7), aparecem como
x
= −αL ϕ̃1/2,j , i = 1, j = 1 : J,
J˜1/2,j
(3.11)
x
= αR ϕ̃I+1/2,j , i = I, j = 1 : J,
J˜I+1/2,j
(3.12)
y
Jˆi,1/2
= −αB ϕ̂i,1/2 , i = 1 : I, j = 1,
(3.13)
y
Jˆi,J+1/2
= αT ϕ̂i,J+1/2 , i = 1 : I, j = J.
(3.14)
3.2 Equações de Balanço espacial
Obtemos as equações de balanço espacial de ordem zero aplicando o operador
32
∫
∫
xi+1/2 yj+1/2
1
hxi hyj
•dxdy ,
(3.15)
xi−1/2 yj−1/2
no sistema de equações de difusão para o nodo
Γij , ou seja nas eq.
(3.8, 3.9 e 3.10), o resultado
aparece na forma
1 ˆy
1
1 ˜x
y
x
(Ji+1/2,j − J˜i−1/2,j
)+
(Ji,j+1/2 − Jˆi,1/2−j
) + Σaij ϕ̄ij =
νΣfij ϕ̄ij ,
hx i
hyj
kef
(3.16)
Dij
J¯ijx = −
(ϕ̃i+1/2,j − ϕ̃i−1/2,j ) ,
hxi
(3.17)
Dij
(ϕ̂i,j+1/2 − ϕ̂i,j−1/2 ) ,
J¯ijy = −
hyj
(3.18)
ϕ̄ij ,
onde denimos o uxo médio no interior do nodo
nodo para cada uma das direções espaciais,
∫
J¯ijx
e
J¯ijy ,
e as correntes médias no interior do
segundo as seguintes expressões
∫
xi+1/2 yj+1/2
1
ϕ̄ij =
hxi hyj
ϕ(x, y)dxdy ,
(3.19)
J x (x, y)dxdy ,
(3.20)
xi−1/2 yj−1/2
∫
∫
xi+1/2 yj+1/2
1
J¯ijx =
hxi hyj
xi−1/2 yj−1/2
∫
∫
xi+1/2 yj+1/2
J¯ijy
1
=
hxi hyj
J y (x, y)dxdy ,
(3.21)
xi−1/2 yj−1/2
também denimos as grandezas médias nas arestas verticais do nodo, aresta esquerda e direita,
na forma
∫
yj+1/2
R̃i±1/2,j
1
=
hyj
Ri±1/2 (y)dy ,
onde, R
= J x ou ϕ ,
(3.22)
yj−1/2
ademais denimos as grandezas médias nas arestas horizontais do nodo, aresta inferior e
superior, na forma
∫
xi+1/2
R̂i,j±1/2
1
=
hxi
Rj±1/2 (y)dy ,
onde, R
= J y ou ϕ ,
(3.23)
xi−1/2
uma representação para um nodo arbitrário das grandezas médias no interior do nodo
ϕ̄ij ,
33
J¯ijx , J¯ijy ;
das grandezas médias nas arestas verticais
nas arestas horizontais
y
Jˆi,j±1/2
, ϕ̂i,j±1/2
x
J˜i±1/2,j
, ϕ̃i±1/2,j ,
e das grandezas médias
é mostrada na gura 3.4.
Figura 3.4: Representação das grandezas associadas a um nodo arbitrário
Γij
.
O sistema de equações de balanço espacial, dado pelas eqs. (3.16, 3.17 e 3.18), fornece
um sistema de três equações lineares e algébricas para cada nodo. Considerando uma grade
espacial de
I×J
nodos, as equações de balanço espacial de cada nodo, e as condições de
contorno do problema teremos o seguinte balanço de incógnitas e equações
Incógnitas médias no nodo (ϕ̄ij ,
J¯ijx , J¯ijy )
x
Incógnitas nas arestas verticais (J˜i±1/2,j ,
3IJ
ϕ̃i±1/2,j )
y
Incógnitas nas arestas horizontais (Jˆi,j±1/2 ,
ϕ̂i,j±1/2 )
2(I + 1)J
2I(J + 1)
Total de incógnitas
7IJ + 2(I + J)
Equações de balanço espacial (3.16, 3.17 e 3.18)
3IJ
Equações dos contornos esquerdo e direito (3.11 e 3.12)
2J
Equações dos contornos superior e inferior (3.13 e 3.14)
2I
Total de equações
3IJ + 2(I + J)
Para auxiliar na visualização das incógnitas presentes na grade espacial ilustramos na gura
3.5 uma discretização em uma malha de
3 × 3.
O balanço anterior nos mostra que o conjunto
de equações de balanço e as correspondentes condições de contorno não fornecem solução
única para o problema de autovalor de difusão. Uma vez que o número de incógnitas supera
o número de equações, teremos um sistema solúvel e indeterminado. Para garantir a solução
34
Figura 3.5: Grade espacial e incógnitas, considerando uma malha de
3 × 3; (•)
grandezas médias
na célula espacial 3 por nodo; (◦) grandezas em arestas verticais 2 por aresta; e (×) grandezas em
arestas horizontais 2 por aresta.
única do problema, devemos introduzir
4IJ
equações auxiliares; ou seja, 4 equações auxiliares
para cada nodo. As equações auxiliares introduzidas, denem as aproximações e o tipo de
método numérico que será utilizado. Nas próximas seções apresentamos alguns dos métodos
numéricos desenvolvidos para problemas de autovalor de difusão.
3.3 Método Diamond Dierence (DD)
O método Diadmond Dierence (DD) foi desenvolvido para problemas de transporte
na formulação de ordenadas discretas [Lewis e Miller 1993], posteriormente ele foi aplicado a
outras formulações de transporte. As equações auxiliares do método DD calculam as grandezas médias no nodo como a média das grandezas nas arestas do nodo. Para o problema
bidimensional de difusão as equações auxiliares do método DD aparecem na forma
ϕ̄ij =
ϕ̃i+1/2,j + ϕ̃i−1/2,j
,
2
(3.24)
ϕ̄ij =
ϕ̂i,j+1/2 + ϕ̂i,j−1/2
,
2
(3.25)
J¯ijx =
J¯ijy
=
x
x
J˜i+1/2,j
+ J˜i−1/2,j
2
y
y
Jˆi,j+1/2
+ Jˆi,j−1/2
2
,
(3.26)
.
(3.27)
35
As aproximações introduzidas, pelo método DD, supõem que as grandezas uxo e
corrente de nêutrons apresentam um comportamento linear no interior do nodo. Esta aproximação impõe restrições ao tamanho dos nodos, que podem ser utilizados na discretização
espacial, à medida que utilizarmos malhas mais nas, a aproximação é mais consistente. A
inuência do tamanho do nodo na qualidade da aproximação pode ser vericada na gura 3.6.
Por este motivo o método DD é classicado como um método de malha na.
Figura 3.6: Aproximação do método DD para o uxo e a corrente no interior do nodo, a esquerda
discretização de malha grossa, a direita discretização de malha na.
As equações de balanço espacial (3.16, 3.17 e 3.18), em conjunto com as equações auxiliares (3.24, 3.25, 3.26 e 3.27) e as condições de contorno (3.11, 3.12, 3.13 e 3.14) fornecem um
sistema de equações solúvel e determinado para o problema de autovalores de difusão. Entretanto, dado o elevado número de equações e incógnitas do sistema é recomendável realizar
transformações algébricas no sistema que nos permitam simplicá-lo.
Começamos substi-
tuindo as equações auxiliares nas equações de balanço. Primeiro, substituímos a eq. (3.24)
no membro esquerdo da equação (3.16). O resultado aparece na forma
1 ˜x
1 ˆy
y
x
(Ji+1/2,j − J˜i−1/2,j
)+
(J
− Jˆi,j−1/2
)+
hx i
hyj i,j+1/2
Σaij
1 νΣfij
ϕ̄ij .
(ϕ̃i+1/2,j + ϕ̃i−1/2,j ) =
2
kef 2
(3.28)
Em seguida, substituímos a eq. (3.25) no membro esquerdo da equação (3.16) e obtemos
1 ˆy
1 ˜x
y
x
(Ji+1/2,j − J˜i−1/2,j
)+
(Ji,j+1/2 − Jˆi,j−1/2
)+
hx i
hyj
1 νΣfij
Σaij
(ϕ̂i,j+1/2 + ϕ̂i,j−1/2 ) =
ϕ̄ij .
2
kef 2
(3.29)
Ao substituirmos a equação (3.26) na equação (3.17) obtemos
2Dij
x
x
=−
+ J˜i−1/2,j
J˜i+1/2,j
(ϕ̃i+1/2,j − ϕ̃i−1/2,j ) .
hx i
(3.30)
36
Por último, substituímos a eq. (3.27) na eq. (3.18) e o resultado aparece na forma
2Dij
y
y
Jˆi,j+1/2
+ Jˆi,j−1/2
=−
(ϕ̂i,j+1/2 − ϕ̂i,j−1/2 ) .
hy j
(3.31)
A varredura numérica para resolver o problema compreende um esquema explícito de
direções alternadas.
x
(J˜i±1/2,j ,
ϕ̃i±1/2,j )
Neste esquema agrupamos as grandezas médias nas arestas verticais
de cada la horizontal conforme gura 3.7. O sistema de equações corres-
pondente a cada linha
j
inclui uma condição de contorno à esquerda eq. (3.11);
2I
equações
modicadas de balanço dadas pelas eqs. (3.28 e 3.30), e uma equação do contorno à direita
(3.12). Portanto, escrevemos
J˜1/2,j + αL ϕ1/2,j = 0 ,
(3.32)
1 ˜
Σaij
(ϕ̃i+1/2,j + ϕ̃i−1/2,j ) =
(Ji+1/2,j − J˜i−1/2,j ) +
hxi
2
1
1 ˆy
y
νΣfij (ϕ̄ij ) −
(J
− Jˆi,j−1/2
),
kef
hyj i,j+1/2
(3.33)
2Dij
x
x
(J˜i+1/2,j
+ J˜i−1/2,j
)+
(ϕ̃i+1/2,j − ϕ̃i−1/2,j ) = 0 ,
hx i
(3.34)
x
J˜I+1/2
− αR ϕ̃I+1/2,j = 0 .
(3.35)
Aqui, no lado esquerdo do sistema aparecem as incógnitas de uxo e corrente nas arestas
verticais e no lado direito o termo de fonte de ssão e o termo de fuga transversal na direção
y.
O termo de fuga transversal é computado com as correntes médias nas arestas horizontais
que serão calculadas com o sistema de equações (3.363.39).
Figura 3.7: Representação das incógnitas do problema nas arestas verticais para uma linha
Por outro lado, os sistemas correspondentes a cada la
x
uxo e corrente médios nas arestas horizontais (Jˆi±1/2,j ,
ϕ̂i±1/2,j ).
i
j.
agrupam as incógnitas de
Uma representação de uma
la vertical e as correspondentes incógnitas é oferecida na gura 3.8. O sistema de equações
correspondente a cada coluna
i
inclui uma condição de contorno inferior equações (3.13),
2I
equações de balanço equações (3.29 e 3.31), e uma equação do contorno superior (3.14). O
37
resultado aparece na forma
y
Jˆi,1/2
= −αB ϕ̂i,1/2 ⇒ Jˆi,1/2 + αB ϕ̂i,1/2 = 0 ,
1 ˆ
Σaij
(Ji,j+1/2 − Jˆi,j−1/2 ) −
(ϕ̂i,j+1/2 + ϕ̂i,j−1/2 ) =
hyj
2
1
1 ˜x
x
)
− J˜i−1/2,j
νΣfij (ϕ̄ij ) −
(J
kef
hxi i+1/2,j
(3.36)
(3.37)
2Dij
y
y
(Jˆi,j+1/2
+ Jˆi,j−1/2
)+
(ϕ̂i,j+1/2 − ϕ̂i,j−1/2 ) = 0
hyj
(3.38)
y
JˆJ+1/2
− αT ϕ̂i,J+1/2 = 0 .
(3.39)
Neste sistema de equações as incógnitas do problema (uxo e corrente) aparecem do lado
esquerdo, e os termos do lado direito dependem da fonte de ssão ou das fugas transversais
na direção
x.
Estes termos de fuga transversal são calculados com as correntes obtidas no
sistema de equações (3.323.35).
Figura 3.8: Representação das incógnitas do problema nas arestas horizontais para uma linha
i.
O algoritmo para resolver o problema de autovalor de difusão inclui iterações externas
e internas. Nas iterações externas utilizamos o tradicional método da potência, para calcular o
autovalor dominante. Nas iterações internas usamos um método numérico para resolver SELA
e um esquema de direções alternadas para determinar os uxos e as correntes nas arestas dos
nodos. A cada iteração interna os valores dos termos de fuga transversal são atualizados.
Os passos detalhados do algoritmo utilizado na implementação do método DD são
38
oferecidos a seguir:
1. Entrada de dados e condições de contorno:D ,
2. Cálculo de auxiliares:
3. Estimativa inicial:
Σa, νΣf , Lx , Ly , I , J , αL,B,R,T , ξϕ , ξke f .
hxi , hyj .
J˜i±1/2,j , Jˆi,j±1/2 , ϕ̃i±1/2,j , ϕ̂i,j±1/2
e
kef .
4. Cálculo da integral da fonte:
- Calcular
ϕ̄ij
usando a equação auxiliar ;
- Calcular a integral da fonte:
IF0 =
∑ ∑
i
j
5. Início do processo iterativo externo (IterExt
Fy =
para
1 ˆy
y
(Ji,j+1/2 − Jˆi,j−1/2
)
hxi
x
x
J˜i±1/2
, ϕ̃i±1/2 ,
usando as equações (3.32), (3.33), (3.34), (3.35)
j = 1 : J.
8. Calcular a fuga na direção
x.
Fx =
9. Obter as estimativas para
(3.38), (3.39) para
1 ˜x
x
(J
− J˜i−1/2,j
)
hxi i+1/2,j
y
y
Jˆj±1/2
, ϕ̂j±1/2
usando o sistema de equações (3.36), (3.37),
i = 1 : I.
10. Vericar convergência no uxo
11. Calcular a integral da fonte
(ξϕ ),
ϕ̄ij usando a equação
∑ ∑ ∑
IF1 = i j ν fij ϕ̄ij hxi hyj .
12. Calcular nova estimativa do
se necessário voltar ao passo 6.
(IF1 ):
- Calcular
-
= 0).
y
6. Cálculo das fugas na direção
7. Obter estimativas para
νΣfij ϕ̄ij hxi hyj .
auxiliar.
kef
(l)
(l−1) IF1
kef = kef
.
IF0
39
13. .
Vericar o critérios de convergência das iterações externas
k (l−1) , IF0 = IF1 , IterExt = IterExt + 1
14. Imprimir resultados,
kef ,
se necessário
(l)
kef =
e voltar ao passo 6.
y
x
J˜i±1/2
, ϕ̃i±1/2 , ϕ̂j±1/2 , kef .
, Jˆj±1/2
3.4 Método Espectro Nodal de Difusão com Aproximação Constante (SNDCN)
O método SND-CN utiliza um processo de discretização espacial semelhante ao descrito
na seção 3.1, ademais utiliza as equações de balanço espacial (3.16, 3.17 e 3.18) obtidas na
seção 3.2.
Para obtermos as equações auxiliares do método, seguimos o procedimento de
integração transversal comum a todos os métodos nodais.
O procedimento de integração
transversal consiste em integrar o sistema de equações de difusão no nodo
Γij
para cada uma
das direções espaciais separadamente.
Para obter as equações integradas transversalmente na direção
y,
aplicamos nas eq.
(3.8), (3.9) e (3.10) o operador
1
hy j
∫
yj+1/2
• dy
yj−1/2
e o resultado aparece na forma
d ˜x
1
1 y
y
Jj (x) + Σaij ϕ̃j (x) =
νΣfij ϕ̃j (x) −
[J
(x) − Jj−1/2
(x)] ,
dx
kef
hyj j+1/2
(3.40)
d
J˜jx (x) = −Dij ϕ̃j (x) ,
dx
(3.41)
J˜jy (x) = −Dij [ϕj+1/2 (x) − ϕj−1/2 (x)] ,
(3.42)
onde denimos as grandezas médias na direção
1
ϕ̃j (x) =
hy j
1
J˜jx (x) =
hy j
1
J˜jy (x) =
hyi
∫
y
como
yj+1/2
ϕ(x, y)dy ,
yj−1/2
∫
yj+1/2
J x (x, y)dy ,
yj−1/2
∫
yj+1/2
yj−1/2
J y (x, y)dy .
40
De forma semelhante, para obter as equações integradas transversalmente na direção
x,
aplicamos nas eq. (3.8), (3.9) e (3.10) o operador
1
hx i
∫
xi+1/2
• dx
xi−1/2
o resultado aparece na forma
1 x
d ˆy
1
x
νΣfij ϕ̂i (y) −
Ji (y) + Σaij ϕ̂i (y) =
[J
(y) − Ji+1/2
(y)] ,
dy
kef
hxi i−1/2
(3.43)
−Dij
[ϕi+1/2 (y) − ϕi−1/2 (y)] ,
Jˆiy (y) =
hxi
(3.44)
d
Jˆiy (y) = −Dij ϕ̂i (y) ,
dy
(3.45)
onde as grandezas médias na direção
x
são denidas como
1
ϕ̂i (y) =
hxi
Jˆix (y)
1
=
hxi
1
Jˆiy (y) =
hxi
∫
xi+1/2
ϕ(x, y)dx ,
xi−1/2
∫
xi+1/2
J x (x, y)dx ,
xi−1/2
∫
xi+1/2
J y (x, y)dx .
xi−1/2
O procedimento de integração transversal permite transformar o sistema de equações
de difusão do nodo
Γij
representado por um sistema de equações diferenciais parciais em
dois problemas de equações diferencias ordinárias.
Este procedimento divide o problema
bi-dimensional original em dois problemas unidimensionais.
O problema dependente de
é formado pelas equações (3.40) e (3.41), e o problema dependente de
equações (3.43) e (3.44).
y
x
é formado pelas
Neste ponto, destacamos que os problemas unidimensionais de
EDOs associados a cada uma das direções espaciais encontram-se acoplados pelos termos de
fuga transversal. Os termos de fuga transversal são identicados no último termo à direita
da Eq. (3.40) e no último termo à direita da eq. (3.43) para cada problema.
Na tabela 3.1 ilustramos um balanço do número de equações e do número de incógnitas para os problemas unidimensionais em cada direção. Observamos que os problemas são
indeterminados, ou seja, o número de incógnitas do sistema de EDOs supera o número de
equações. Para podermos resolver estes problemas, aproximações simplicadoras devem ser
41
introduzidas.
Tabela 3.1: Balanço de equações e incógnitas para os problemas unidimensionais de EDOs gerados
pela integração transversal.
Problema x
Equações: (3.40), (3.41)
Problema y
Equações: (3.43), (3.44)
y
x
Incógnitas: J˜i (y), ϕ̃j (x), Ji±1/2 (y)
Incógnitas:
y
J˜j (x), ϕ̃j (x), Jj±1/2
(x)
Até o momento, as equações integradas transversalmente representam um problema
exato, neste ponto introduzimos as aproximações do método SND-CN, que consistem em
aproximar por polinômios de ordem zero os termos de fuga transversal. O resultado aparece
na forma
y
Jj±1/2
(x)
1
∼
=
hx i
x
Ji±1/2
(y)
1
∼
=
hyj
∫
xi+1/2
xi−1/2
∫
yj+1/2
y
y
Jj±1/2
(x)dx ∼
= Jˆi,j±1/2 ,
x
x
Ji±1/2
(y)dy ∼
,
= J˜i±1/2,j
yj−1/2
a aproximação constante nos termos de fuga transversal, estabelece que as correntes nêutronicas nas arestas dos nodos são aproximadas pelos valores médios nessas arestas. Reescrevemos
então os últimos termo à direita das eq. (3.40), (3.43) na forma
[
y
y
Jj+1/2
(x) − Jj−1/2
(x)
]
hyj
[
x
x
Ji+1/2
(y) − Ji−1/2
(y)
hxi
]
1 ˆy
y
y
∼
(J
− Jˆi,j−1/2
)∼
=
= F̂i,j ,
hyj i,j+1/2
(3.46)
1 ˜x
x
x
∼
(Ji+1/2,j − J˜i−1/2,j
)∼
.
=
= F̃i,j
hy j
(3.47)
Substituindo as aproximações (3.46) e (3.47) nas eq.
(3.40) e (3.43) obtemos as equações
unidimensionais integradas transversalmente do método espectronodal de difusão com aproximação constante, cujo o resultado aparece na forma
1
d ˜x
y
Jj (x) + Σaij ϕ̃j (x) =
νΣfij ϕ̃j (x) − F̂i,j
,
dx
kef
(3.48)
d
J˜jx (x) = −Dij ϕ̃j (x) ,
dx
(3.49)
1
d ˜y
x
,
Ji (y) + Σaij ϕ̃i (y) =
νΣfij ϕ̂i (y) − F̃i,j
dy
kef
(3.50)
42
d
Jˆiy (y) = −Dij ϕ̂i (y) .
dy
(3.51)
As eq. (3.48) e (3.49) representam o sistema de equações integradas transversalmente
do método SND-CN para o problema dependente de
x;
e as eq. (3.50) e (3.51) representam
o sistema de equações integradas transversalmente do método SND-CN para o problema
dependente de
y.
Para cada um destes problemas por separado pode ser encontrada uma
solução analítica.
Em seguida resolvemos analiticamente o problema de EDOs representado pelas eq.
(3.48) e (3.49). Começamos, substituindo a eq. (3.49) na eq. (3.48) e obtemos
−Dij
1
d2 ϕ̃j (x)
y
+ Σai j ϕ̃j (x) =
νΣfij ϕ̃j (x) − F̂i,j
,
2
dx
kef
(3.52)
a expressão (3.52) representa uma EDO não homogênea com coecientes constantes, e sua
solução pode ser calculada como uma combinação linear da solução homogênea e uma solução
particular que independe da variável
x,
ou seja
ϕ̃j (x) = ϕ̃hj (x) + ϕ̃pj .
Aqui os sobrescritos
h
e
p
referem-se as componentes homogênea e particular da solução
respectivamente. Substituindo a componente particular da solução na eq. (3.52) e agrupando
apropriadamente obtemos
ϕ̃pj = (
F̂ijy
1
νΣfij − Σaij
kef
).
Para calcularmos a componente homogênea da solução propomos o
(3.53)
ansatz
x
ϕ̃hyj (x) = e ξ ,
em seguida substituímos o
ansatz
na eq. (3.52) e agrupando apropriadamente obtemos os
autovalores do problema na forma
v
u
u
ξℓ = u
t
Dij
, com ℓ = 1, 2 ,
1
Σaij −
νΣfij
kef
a componente homogênea da solução aparece na forma
(3.54)
43
ϕ̃hj (x)
=
2
∑
( ξχ )
αℓ e
, ℓ = 1, 2 .
ℓ
(3.55)
ℓ=1
Assim obtemos a solução geral analítica da eq.
transversalmente na direção
y
(3.52), que corresponde ao uxo integrado
no interior do nodo
ϕ̃j (x).
Somando a solução homogênea eq.
(3.55) e a solução particular eq. (3.53), o resultado aparece na forma
ϕ̃j (x) =
2
∑
( ξχ )
αℓ e
ℓ
+(
ℓ=1
onde
αℓ
são constantes arbitrárias e
ξℓ
F̂ijy
1
νΣfij − Σaij
kef
) , ℓ = 1, 2 ,
(3.56)
são os autovalores do problema (3.52) representados
pela equação (3.54).
Obtemos a solução geral analítica para a corrente neutrônica média na direção
[J˜jx (x)],
substituindo a eq.
(3.56) na eq.
(3.49).
y
Calculando a primeira derivada o resul-
tado aparece como
J˜jx (x) = −Dij
2
∑
αℓ
ℓ=1
ξℓ
χ
e ξℓ .
(3.57)
As eq. (3.56) e (3.57) representam as soluções gerais analíticas do sistema de equações
integradas transversalmente dependente de
x.
Seguindo um procedimento semelhante podem
ser obtidas as soluções gerais analíticas do problema dependente de
ϕ̂i (y)
e
Jˆiy (y)
2
∑
( ξy )
βℓ e
ℓ
+(
ℓ=1
Jˆiy (y) = −Dij
2
∑
βℓ
ℓ=1
βℓ
ou seja determinar
a partir das eq. (3.50) e (3.51). Neste caso as soluções analíticas são:
ϕ̂i (y) =
onde
y,
representa constantes arbitrárias e
ξℓ
ξℓ
F̃ x
1
νΣfi j − Σaij
kef
e
( ξy )
ℓ
, ℓ = 1, 2 ,
),
(3.58)
(3.59)
os autovalores do problema dependente de
que coincidem com os autovalores do problema dependente de
y,
x e são calculados usando a eq
(3.54).
Neste ponto devemos construir um esquema numérico convergente que preserve incondicionalmente as soluções gerais analíticas das equações integradas transversalmente do
método SND-CN. Ou seja, nosso esquema numérico deverá preservar as soluções dadas pelas
eq. (3.56) a (3.59) no interior de cada nodo. Ademais o esquema numérico deve incluir: as
condições de contorno do problema, e condições de continuidade do uxo e da corrente nas
44
interface entre os nodos.
Conforme foi discutido na seção 3.2, as equações de balanço em conjunto com as condições de contorno fornecem um sistema de equações indeterminado, para garantir a unicidade
da solução, 4 equações auxiliares para cada nodo devem ser adicionadas ao problema.
No
método SND-CN as equações auxiliares calculam as grandezas médias no interior do nodo
através de médias ponderadas das grandezas nas arestas do nodo.
Portanto, as equações
auxiliares do método SND-CN são escritas como
ϕ̄i,j =
γij
(ϕ̃i+1/2,j − ϕ̃i−1/2,j ) + χ̃ij ,
2
(3.60)
ϕ̄i,j =
λij
(ϕ̂i,j+1/2 − ϕ̂i,j−1/2 ) + χ̂ij ,
2
(3.61)
γij ˜x
x
J¯ijx =
(Ji+1/2,j + J˜i−1/2,j
),
2
(3.62)
λij ˆy
y
J¯ijy =
(Ji,j+1/2 + Jˆi,j−1/2
).
2
(3.63)
γi,j , λi,j , χ̃i,j
Aqui os valores das constantes
e
χ̂i,j
devem ser calculados de forma que as
soluções gerais analíticas das equações integradas transversalmente eqs. (3.56) a (3.59) sejam
preservadas incondicionalmente no interior do nodo.
Especicamente a constante
preservar a componente homogênea das soluções do problema dependente de
eq. (3.56) e (3.57). A constante
χ̃i,j
χ̂i,j ,
λi,j
deve
dadas pelas
deve preservar a componente particular da solução dada
pela eq. (3.56). Em relação ao problema dependente de
eqs. (3.58) e (3.59), a constante
x
γi,j
y,
soluções analíticas dadas pelas
deve preservar a componente homogênea e a constante
a componente particular.
Para calcularmos as constantes que aparecem nas equações auxiliares, utilizamos as
soluções analíticas gerais para obter expressões para as grandezas
x
J¯ijy , J˜i±1/2,j
e
y
Jˆi,j±1/2
.
ϕ̄i,j , ϕ̃i±1/2,j , ϕ̂i,j±1/2 , J¯ijx ,
Em seguida, substituímos essas expressões nas equações auxiliares e
após um pouco de álgebra obtemos as seguintes expressões
y
(γi,j − 1)F̂i,j
,
χ̃ij =
hyj ΣTi,j
(3.64)
45
γij =

( )
hx i
2ξ

ℓ

tgh
,



2ξℓ
 hx i
ξℓ ∈ ℜ ,

( )


hx i
2θℓ


tg
,

hx i
2θℓ
√
ξℓ = îθ , θ ∈ ℜ , î = −1 ,
(3.65)
x
(λi,j − 1)F̃i,j
χ̂ij =
,
hxi ΣTi,j
( )

hy j
2ξℓ


,
tgh


2ξℓ

 hy j
λij =

( )


hy j
2θℓ


tg
,

hy j
2θℓ
(3.66)
ξℓ ∈ ℜ ,
(3.67)
ξℓ = îθ , θ ∈ ℜ , î =
√
−1 ,
onde denimos
ΣTij = Σai,j −
1
Σfi,j .
kef
(3.68)
Uma vez determinados os parâmetros das equações auxiliares, devemos construir os
sistemas lineares e algébricos associados ao método SND-CN, e propor um esquema iterativo
que permita resolver os sistemas e obtermos uma aproximação numérica da solução do problema. Para isso, devemos substituir as equações auxiliares, eq. (3.60) a (3.63), nas equações
do balanço espacial, eq. (3.16) a (3.18), o que nos permite eliminar como incógnitas as grandezas médias nos nodos (ϕ̄i,j ,
J¯ijx
e
J¯ijy ) e obtermos sistemas de equações que envolvam apenas
as grandezas médias nas arestas dos nodos.
Começamos substituindo a equação auxiliar (3.60) na equação de balanço (3.16) o
resultado aparece como
x
x
J˜i+1/2,j
− J˜i−1/2,j
hx i
+ γij
y
y
(Jˆi,j+1/2
− Jˆi,j−1/2
)
hyj
+
) νΣf γ (
)
Σaij γij (
ij ij
ϕ̃i+1/2,j + ϕ̃i−1/2,j =
ϕ̃i+1/2,j + ϕ̃i−1/2,j .
2
2kef
(3.69)
Em seguida, substituímos a equação auxiliar (3.61) na equação de balanço (3.16) e obtemos
46
λij
x
x
(J˜i+1/2,j
− J˜i−1/2,j
)
hx i
)
+
y
y
(Jˆi,j+1/2
− Jˆi,j−1/2
)
hy j
+
)
Σaij λij (
νΣfij λij (
ϕ̂i,j+1/2 + ϕ̂i,j−1/2 =
ϕ̂i,j+1/2 + ϕ̂i,j−1/2 ;
2
2kef
(3.70)
além disso, substituindo as equações auxiliares (3.62) e (3.63) nas equações de balanço (3.17)
e (3.18) respectivamente obtemos,
(
(
As eq.
)
)
2Dij (
x
x
J˜i+1/2,j
+ J˜i−1/2,j
=−
ϕ̃i+1/2,j − ϕ̃i−1/2,j
γij hxi
(3.71)
)
)
2Dij (
y
y
Jˆi,j+1/2
+ Jˆi,j−1/2
=−
ϕ̂i,j+1/2 − ϕ̂i,j−1/2 .
λij hyj
(3.72)
(3.69, 3.70, 3.71 e 3.72) representam um sistema de equações onde apenas
grandezas médias nas arestas dos nodos estão envolvidas.
Ainda podemos simplicar este,
sistema de equações eliminando as correntes nêutronicas, para isso, escrevemos o sistema
em notação matricial e resolvemos em relação às correntes utilizando a regra de Cramer. O
resultado para as respectivas correntes em cada aresta do nodo aparece na forma
[
x
J˜i+1/2,j
]
Dij
=
−
ϕ̃i+1/2,j +
4(1 − γij λij )
hxi γij
]
[
−hxi γij (Σaij − k1ef νΣf )
Dij
+
ϕ̃i−1/2,j +
4(1 − γij λij )
hxi γij
]
[
)
hxi γij λij (Σaij − k1ef νΣf ) (
ϕ̂i,j+1/2 + ϕ̂i,j−1/2
4(1 − γij λij )
−hxi γij (Σaij −
[
1
νΣf )
kef
]
D
ij
x
=
J˜i−1/2,j
−
ϕ̃i+1/2,j +
4(1 − γij λij )
hxi γij
]
[
hxi γij (Σaij − k1ef νΣf )
Dij
ϕ̃i−1/2,j +
+
4(1 − γij λij )
hxi γij
[
]
)
hxi γij λij (Σaij − k1ef νΣf ) (
ϕ̂i,j+1/2 + ϕ̂i,j−1/2
4(1 − γij λij )
hxi γij (Σaij −
(3.73)
1
νΣf )
kef
(3.74)
47
[
−hyj λij (Σaij −
1
νΣf )
kef
]
Dij
ϕ̂i,j+1/2 +
4(1 − γij λij )
hyj λij
[
]
hyj λij (Σaij − k1ef νΣf )
Dij
−
ϕ̂i,j−1/2 +
4(1 − γij λij
hyj λij
]
[
)
hyj γij λij (Σaij − k1ef νΣf ) (
ϕ̃i+1/2,j + ϕ̃i−1/2,j
4(1 − γij λij )
y
Jˆi,j+1/2
=
[
−
(3.75)
]
D
ij
y
Jˆi,j−1/2
=
ϕ̂i,j+1/2 +
−
4(1 − γij λij )
hyj λij
[
]
hyj λij (Σaij − k1ef νΣf )
Dij
+
ϕ̂i,j−1/2 +
4(1 − γij λij
hyj λij
[
]
)
−hyj γij λij (Σaij − k1ef νΣf ) (
ϕ̃i+1/2,j + ϕ̃i−1/2,j
4(1 − γij λij )
hyj λij (Σaij −
1
νΣf )
kef
(3.76)
Em seguida podemos realizar um procedimento algébrico com as equações resultantes
x
que permita eliminar as correntes nas arestas dos nodos (J˜i±1/2,j e
y
Jˆi,j+1/2
);
obtendo assim
um sistema onde as incógnitas são apenas os uxos nas arestas dos nodos (ϕ̃i±1/2,j e
ϕ̂i,j±1/2 ).
Construímos então um sistema de equações para as arestas verticais onde as incógnitas
são apenas os uxos (ϕ̃i±1/2,j ). As equações (3.73 e 3.74) representam a corrente nas arestas
verticais internas dos nodos e as eq. (3.11 e 3.12) representam as correntes nas arestas verticais
das fronteiras esquerda e direita respectivamente. Aplicamos a continuidade da corrente nas
arestas verticais internas, ou seja, fazemos
eq. (3.73). O resultado aparece na forma
i = i+1
na eq. (3.74) e igualamos o resultado à
48
[
]
−Dij
hxi γij Σaij
+
ϕ̃i−1/2,j +
hxi γij 4(1 − γij λij )
[
]
hxi+1 γi+1,j Σai+1,j
Dij
Di+1,j
+
+
ϕ̃i+1/2,j +
hxi γij hxi+1 γi+1,j 4(1 − γi+1,j λi+1,j )
[
]
hxi+1 γi+1,j Σai+1,j
Di+1,j
−
+
ϕ̃i+3/2,j ≡
hxi+1 γi+1,j 4(1 − γi+1,j λi+1,j )
[
]
hxi γij k1ef νΣfij
ϕ̃i−1/2,j +
4(1 − γij λij )
]
[
hxi γij k1ef νΣfij hxi+1 γi+1,j k1ef νΣfi+1,j
+
ϕ̃i+1/2,j +
4(1 − γij λij )
4(1 − γi+1,j λi+1,j )
]
[
hxi+1 γi+1,j k1ef νΣfi+1,j
ϕ̃i+3/2,j +
4(1 − γi+1,j λi+1,j )
)
hxi γij λij ΣTij (
ϕ̂i,j+1/2 + ϕ̂i,j−1/2 +
4(1 − γij λij )
)
hxi+1 γi+1,j λi+1,j ΣTi+1,j (
ϕ̂i,j+1/2 + ϕ̂i,j−1/2
4(1 − γi+1,j λi+1,j )
(3.77)
Em seguida consideramos a condição da continuidade da corrente no contorno esquerdo, ou
seja consideramos
i = 1
na eq.
(3.74) e igualamos o resultado à eq.
(3.11).
Agrupando
apropriadamente obtemos
[
]
D1,j
h1 γ1,j Σa1,j
L
+
+ α ϕ̃1/2,j +
h1,j γ1,j 4(1 − γ1,j λ1,j )
]
[
−D1,j
h1 γ1,j Σa1,j
ϕ̃3/2,j ≡
+
h1,j γ1,j 4(1 − γ1,j λ1,j )
[
]
h1 γ1,j k1ef νΣf
ϕ̃3/2,j +
4(1 − γ1,j λ1,j )
[
]
h1 γ1,j k1ef νΣf
ϕ̃1/2,j +
4(1 − γ1,j λ1,j )
[
]
)
h1 γ1,j λ1,j ΣT1,j (
ϕ̂1,j+1/2 + ϕ̂1,j−1/2
4(1 − γ1,j λ1,j )
(3.78)
Finalmente consideramos a condição da continuidade da corrente no contorno direito, ou seja
fazemos
i=I
na eq. (3.73) e igualamos o resultado a eq. (3.12), o resultado aparece na forma
49
[
]
−DI,j
hI γI,j ΣaI,j
+
ϕ̃I−1/2,j +
hI,j γI,j 4(1 − γI,j λI,j )
[
]
DI,j
hI γI,j ΣaI,j
R
+
+ α ϕ̃I+1/2,j ≡
hI,j γI,j 4(1 − γI,j λI,j )
[
]
hI γI,j k1ef νΣf
ϕ̃I+1/2,j +
4(1 − γI,j λI,j )
[
]
hI γI,j k1ef νΣf
ϕ̃I−1/2,j +
4(1 − γI,j λI,j )
[
]
)
hI γI,j λI,j ΣTI,j (
ϕ̂I,j+1/2 + ϕ̂I,j−1/2
4(1 − γI,j λI,j )
(3.79)
O sistema de equações formado por (3.77, 3.78 e 3.79) representa um esquema de três
pontos para cada uma das camadas horizontais
j = 1...J
do domínio.
A resolução dos
J
sistemas nos permite calcular todos os uxos (ϕ̃i±1/2,j ) no domínio de cálculo.
Analogamente podemos obter sistemas tridiagonais para cada camada vertical
onde as incógnitas são os uxos (ϕ̂i,j±1/2 ).
i = 1..I
Consideramos as equações (3.75 e 3.76), que
representam as correntes nas arestas horizontais internas dos nodos, além das equações (3.13 e
3.14) que representam as correntes nas arestas horizontais das fronteiras superiores e inferiores,
respectivamente. Deste modo, aplicamos a continuidade da corrente nas arestas horizontais
internas, ou seja, fazemos
obtemos
j = j + 1 na eq.
(3.76) e igualamos o resultado à eq. (3.75). Assim,
50
[
[
Dij
Di,j+1
+
hyj λij hyj+1 λi,j+1
]
hyj λij Σaij
−Dij
+
ϕ̂i,j−1/2 +
hyj λij 4(1 − γij λij )
]
hj+1 λi,j+1 Σai,j+1
hj λi,j Σai,j
+
ϕ̂i,j+1/2 +
+
4(1 − γi,j λi,j ) 4(1 − γi,j+1 λi,j+1 )
[
]
−Di,j+1
hj+1 λi,j+1 Σai,j+1
ϕ̂i,j+3/2 ≡
+
hyj+1 λi,j+1 4(1 − γi,j+1 λi,j+1 )
]
[
hyj λij k1ef νΣfij
ϕ̂i,j−1/2 +
4(1 − γij λij )
[
]
hyj λij k1ef νΣf
hyj+1 λi,j+1 k1ef νΣf
+
ϕ̂i,j+1/2 +
4(1 − γij λij )
4(1 − γi,j+1 λi,j+1 )
[
]
hyj+1 λi,j+1 k1ef νΣf
ϕ̂i,j+3/2 +
4(1 − γi,j+1 λi,j+1 )
)
hyj γij λij ΣTij (
ϕ̃i+1/2,j + ϕ̃i−1/2,j +
4(1 − γij λij )
)
hyj+1 γi,j+1 λi,j+1 ΣTi,j+1 (
ϕ̃i+1/2,j + ϕ̃i−1/2,j .
4(1 − γi,j+1 λi,j+1 )
(3.80)
Em seguida, consideramos a condição de continuidade da corrente no contorno inferior, ou
seja consideramos
j = 1
na eq.
(3.76) e igualamos o resultado à eq.
(3.13).
Agrupando
apropriadamente obtemos
[
]
Di,1
h1 λi,1 Σai,1
B
+
+ α ϕ̂i,1/2 +
h1 λi,1 4(1 − λi,1 γi,1 )
[
]
−Di,1
h1 λi,1 Σai,1
ϕ̂i,3/2 ≡
+
h1 λi,1 4(1 − λi,1 γi,1 )
[
]
h1 λi,1 k1ef νΣi,1
ϕ̂i,1/2 +
4(1 − λi,1 γi,1 )
[
]
h1 λi,1 k1ef νΣfi,1
ϕ̂i,3/2 +
4(1 − λi,1 γi,1 )
[
]
)
h1 λi,1 γi,1 ΣTi,1 (
ϕ̃i+1/2,1 + ϕ̃i−1/2,1
4(1 − λi,1 γi,1 )
(3.81)
Finalmente, consideramos a condição da continuidade da corrente no contorno superior, ou
seja fazemos
forma
j=J
na eq. (3.75) e igualamos o resultado a eq. (3.14), o resultado aparece na
51
[
]
−Di,J
hJ λi,J Σai,J
+
ϕ̂i,J−1/2 +
hJ λi,J
4(1 − λi,J γi,J )
[
]
Di,J
hJ λi,J Σai,J
T
+
+ α ϕ̂i,J+1/2 ≡
hJ λi,J
4(1 − λi,J γi,J )
[
]
hJ λi,J k1ef νΣfi,J
ϕ̂i,J+1/2 +
4(1 − λi,J γi,J )
[
]
hJ λi,J k1ef νΣfi,J
ϕ̂i,J−1/2 +
4(1 − λi,J γi,J )
[
]
)
hJ λi,J γi,J ΣTi,J (
ϕ̃i+1/2,J + ϕ̃i−1/2,J
4(1 − λi,J γi,J )
(3.82)
O sistema de equações formado por (3.80, 3.81 e 3.82) representa um esquema de três
pontos para cada uma das camadas verticais
i = 1...I
do domínio. A resolução dos
I
sistemas
nos permite calcular todos os uxos (ϕ̂i±1/2,j ) no domínio de cálculo.
Os sistemas de equações construídos constituem as equações em diferenças do método
SND-CN. Estes sistemas podem ser combinados em uma varredura que utilize um esquema
explícito de direções alternadas que convirja a solução do problema de autovalor.
Maiores
informações sobre o algoritmo de cálculo e o esquema de direções alternadas do SND-CN
podem ser encontradas em [Dominguez 2000].
4 MÉTODO ESPECTRONODAL DE GRADES COMPOSTAS
Os
trabalhos
[Dominguez e Barros 2007],
[Dominguez e Barros 2010]
e
[Barros et al. 2003] mostraram que os métodos espectronodais são ecientes e precisos para
cálculos de malhas grossas em domínios bidimensionais. Entretanto, as principais desvantagens desta família de métodos são: o alto custo algébrico para obtermos as equações discretizadas de varredura, e o alto custo algorítmico da implementação computacional. Visando
amenizar estas desvantagens e ao mesmo tempo aproveitar a precisão dos métodos espectronodais surgiram os métodos espectronodais de grades compostas. O diferencial dos métodos
de grades compostas é introduzir aproximações para os termos de fuga transversal em nível
de região, conseguindo assim desacoplar as direções espaciais e obter uma solução aproximada
do problema através da solução de dois problemas unidimensionais.
A aproximação de fugas nos métodos de grades compostas é mais grosseira que nos métodos nodais convencionais. Entretanto, o método espectronodal de grades compostas provou
ser eciente para calcular grandezas integrais em problemas de fonte xa e aproximação de
ordenadas discretas [Dominguez et al. 2009]. O principal objetivo desta pesquisa é a construção do método espectronodal de grades compostas, para resolver numericamente problemas de
autovalor de difusão em geometria
X, Y
e aproximação de uma velocidade com aproximação
constante para os termos de fuga transversal na região (CSG-SND-CN, Composite Spatial
Grid - Spectral Nodal Diusion - Constant Nodal).
No método CSG-SND-CN, o processo de discretização das variáveis espaciais é feito
em dois estágios e combina um método não-convencional de grades compostas com o processo
de integração transversal dos métodos nodais de difusão tradicionais [Lawrence 1986].
No
primeiro estágio utilizamos uma grade espacial retangular com células grandes; em domínios
heterogêneos esta grade espacial pode ser coincidente com as regiões materiais dentro do
núcleo do reator nuclear (elementos combustíveis). Dentro desta grade, integramos transversalmente as equações de difusão, separadamente nas direções coordenadas
x e y.
Neste ponto,
introduzimos aproximações constantes para os termos de fuga transversal produzindo dois problemas unidimensionais, acoplados pelos termos de fuga transversal. No segundo estágio de
discretização, usamos uma grade espacial na para discretizar cada problema unidimensional
nodal de difusão, e utilizamos o método SND para resolver os problemas unidimensionais
[Dominguez et al. 2009]. O método SND trata analiticamente os termos de fonte de ssão e
gera soluções numéricas que são completamente livres de erros de truncamento espacial para
53
problemas unidimensionais. Como as direções espaciais são acopladas pelos termos de fuga
transversal, utilizamos uma técnica explícita de direções alternadas para convergir a solução
numérica. Os fundamentos do método CSG-SND-CN são mostrados nas próximas seções.
4.1 Bases Matemáticas do Método CSG-SND-CN
Considerando um domínio retangular homogêneo de altura
Ly
e largura
Lx
como ilus-
trado na gura (4.1). Neste domínio a equação de difusão de nêutrons e aproximação de uma
velocidade é descrita pelas expressões (3.1), (3.2) e (3.3), com as correspondentes condições
de contorno de tipo albedo eq. (3.4-3.7).
Figura 4.1: Domínio da solução e grade espacial grossa exterior
Λ
O primeiro estágio de discretização introduz uma grade espacial externa
por células espaciais largas. Esta grade divide o domínio em
uma destas células tem altura
uma célula espacial
Λij
Hjy
e largura
Hix
I ×J
Λ
células grossas
formada
Λij
cada
como mostrado na gura (4.1). No interior de
as propriedades do meio permanecem constantes, ou seja, as células
da grade espacial externa devem coincidir com as regiões do domínio. A equação de difusão
no interior de cada célula espacial
Λij
é descrita pelas expressões (3.8), (3.9) e (3.10).
Para obter as equações de difusão integradas transversalmente do método CSG-SGFCN sobre a grade
Λ,
aplicamos ao sistema de equações de difusão (3.8), (3.9) e (3.10) o
operador
∫
ut+1/2
1
Htu
onde
u=x
ou
y
e
t=i
ou
j.
• du
ut−1/2
(4.1)
54
Ao considerar
u = y
e
t = j
na eq. (4.1) e aplicando o operador resultante nas eq.
(3.8-3.10), integramos o sistema de equações de difusão na direção
de difusão integradas transversalmente dependentes de
d ˜x
1
Jj (x) + Σaij ϕ̃j (x) =
νΣfij ϕ̃j (x) −
dx
kef
x.
y
e obtemos as equações
O resultado aparece na forma
[
]
y
y
Jj+1/2 (x) − Jj−1/2 (x)
Hjy
,
(4.2)
d
J˜jx (x) = −Dij ϕ̃j (x) ,
dx
(4.3)
]
Dij [
J˜jy (x) = − y ϕj+1/2 (x) − ϕj−1/2 (x) .
Hj
(4.4)
Nas eq. (4.2-4.4) foram denidas as grandezas
1
τ̃j (x) = y
Hj
∫
yj+1/2
τ (x, y)dy
τ = ϕ, J x
para
ou
Jy
e
(x, y) ∈ Λij .
(4.5)
yj−1/2
Se integramos o sistema de equações de difusão na direção
integradas transversalmente dependentes de
y,
x
obtemos as equações
ou seja, consideramos
u = x
e
t = i
no
operador (4.1) e aplicamos o resultado nas eq. (3.8-3.10), obtemos
d ˆy
1
Ji (y) + Σaij ϕ̂i (y) =
νΣfij ϕ̂i (y) −
dy
kef
[
]
x
x
Ji+1/2
(y) − Ji−1/2
(x)
Hix
,
(4.6)
]
Dij [
J˜ix (y) = − x ϕi+1/2 (y) − ϕi−1/2 (y) .
Hi
(4.7)
d
J˜iy (y) = −Dij ϕ̂i (y) ,
dy
(4.8)
Nas eq. (4.6-4.8) foram denidas as grandezas
1
τ̂i (y) = x
Hi
∫
xi+1/2
τ (x, y)dx
para
τ = ϕ, J x
ou
Jy
e
(x, y) ∈ Λij .
(4.9)
xi−1/2
As equações de difusão de malha grossa integradas transversalmente em uma direção
espacial sobre a célula
Λij
formam um sistema de
3
equações diferenciais ordinárias (EDO)
com 4 incógnitas para cada direção espacial. Na direção
x o sistema é formado pelas equações
55
y
J˜jx (x), J˜jy (x), ϕ̃j (x) e Jj±1/2
(x).
(4.2), (4.3) e (4.4) com as incógnitas:
sistema é formado pelas equações (4.6), (4.7) e (4.8) com as incógnitas:
x
Ji±1/2
(y).
y
o
Jˆiy (y), Jˆix (y), ϕ̂i (y)
e
Na direção espacial
Para obtermos solução única nos sistemas de equações integradas transversalmente
para cada direção espacial introduzimos a aproximação de corrente constante nos termos de
fuga transversal ao longo das arestas da célula
Λij .
Ou seja os últimos termos das eq. (4.2) e
(4.6) são aproximados por polinômios de ordem zero, que são escritos na forma
onde
y
Jˆi,j±1/2
y
y
Jj±1/2
(x) ∼
= Jˆi,j±1/2 ,
(4.10)
x
x
Ji±1/2
(y) ∼
,
= J˜i±1/2,j
(4.11)
representa o valor da corrente média na direção
da célula espacial; e
x
J˜i±1/2,j
x
na aresta superior ou inferior
representa o valor da corrente média na direção
y
na aresta
esquerda ou direita da célula espacial.
Substituindo a aproximação (4.10) na Eq. (4.2), obtemos:
d ˜x
1
y
Jj (x) + Σaij ϕ̃j (x) =
νΣfij ϕ̃j (x) − F̂i,j
,
dx
kef
onde os termos de fuga transversal na célula
Λij
ao longo da direção
(4.12)
y
y , F̂i,j
,
foram denidos
como
[
y
F̂i,j
=
y
Jj+1/2
(x)
−
]
y
Jj−1/2
(x)
Hjy
≃
y
y
Jˆi,j+1/2
− Jˆi,j−1/2
Hjy
.
(4.13)
Substituindo a aproximação (4.11) na Eq. (4.6), obtemos:
d ˆy
1
x
,
Ji (y) + Σaij ϕ̂i (y) =
νΣfij ϕ̂i (y) − F̂i,j
dy
kef
onde os termos de fuga transversal na célula
Λij
ao longo da direção
(4.14)
x
,
x, F̂i,j
foram denidos
como
x
F̂i,j
=
[
]
x
x
Ji+1/2 (y) − Ji−1/2 (y)
Hix
.
(4.15)
Neste ponto, temos dois sistema de EDOs de duas equações em duas incógnitas. O
sistema composto pelas Eqs. (4.12) e (4.3) com as incógnitas
espacial
x.
E o sistema formado pelas eq.
J˜jx (x)
e
ϕ̃j (x)
para a direção
(4.14) e (4.7) com as incógnitas
Jˆiy (y)
e
ϕ̂i (y)
56
para a direção espacial
y.
Estes sistemas unidimensionais podem ser resolvidos separada-
mente, no entanto eles estão acoplados pelos termos de fuga transversal. Podemos obter a
solução analítica geral para estes sistemas usando uma técnica de análise espectral, conforme
descrito para o método espectronodal de difusão com aproximação constante (SND-CN) em
[Barros et al. 2003]. As soluções gerais analíticas para o problema dependente de
x aparecem
como
2
∑
ϕ̃j (x) =
(x/ξl )
βl e
l=1
J˜jx (x) = −Dij
y
F̂i,j
−
ΣTij
2
∑
βl
l=1
ξl
(4.16)
e(x/ξ) .
(4.17)
As soluções gerais analíticas para o problema dependente de
ϕ̂i (y) =
2
∑
βl e(x/ξl ) −
l=1
Jˆiy (y) = −Dij
2
∑
βl
l=1
Nas equações acima
βl
ξl
y
aparecem como
x
F̃i,j
,
ΣTij
(4.18)
e(x/ξ) ,
representa constantes arbitrárias,
(4.19)
ξ
são os autovalores do sistema
denidos em [Barros et al. 2004], além disso, denimos a constante
ΣTij = Σaij −
1
νΣfij .
kef
(4.20)
O segundo estágio do processo de discretização consiste em introduzir duas grades
internas nas e independentes
Γx
e
Γy , uma para cada direção espacial, através desta discreti-
zação podemos resolver os problemas unidimensionais separadamente. A grade na
cada célula grossa
Λij
em
M
cada célula grossa
Λij
em
P
nodos de altura
nodos de altura
Hjy
hyi,j
e largura
hxi,j .
A outra grade na
e a correspondente largura
Hix .
Γx
divide
Γy
divide
A geometria
das duas grades espaciais internas é mostrada na gura 4.2.
Neste ponto, obtemos as equações discretizadas de balanço espacial para cada um dos
problemas unidimensionais. Para isso devemos aplicar o operador
∫
ur+1/2
1
hui,j
• du ,
(4.21)
xr−1/2
nos sistemas de equações para cada uma das direções espaciais. Na eq. (4.21)
u=x
ou
y
e
57
Figura 4.2: Grades internas nas na direção
x (Γx ),
e na direção
y (Γy ),
para uma célula grossa
Λij .
r=m
ou
p.
Para obtermos as equações de balanço do problema dependente de
u = x
e
r = m
na eq.
(4.21), aplicamos o operador resultante nas eq.
x
consideramos
(4.12 e 4.3) e o
resultado aparece na forma
1 ˜x
1
y
x
(Jm+1/2,j − J˜m−1/2,j
) + Σaij ϕ̄xm,j =
νΣfij ϕ̄xm,j − F̂i,j
,
x
hi,j
kef
(4.22)
Dij
x
J¯m,j
= − x (ϕ̃m+1/2,j − ϕ̃m−1/2,j ) ,
hi,j
(4.23)
onde foram denidas as grandezas
x
τ̄m,j
para
1
= x
hi,j
xm−1/2 ≤ x ≤ xm+1/2
problema dependente de
y
e
∫
xm+1/2
τ̃j (x)dx,
para
τj = ϕj
ou
Jjx ;
(4.24)
m−1/2
yj−1/2 ≤ y ≤ yj+1/2 .
consideramos
u=y
e
Para obtermos as equações de balanço do
r=p
na eq. (4.21), aplicamos o operador
resultante nas eq. (4.14 e 4.8), o resultado aparece na forma
1 ˆy
1
y
x
ˆy
νΣfij ϕ̄yi,p − F̃i,j
y (Ji,p+1/2 − Ji,p−1/2 ) + Σaij ϕ̄i,p =
hi,j
kef
(4.25)
Dij
y
J¯i,p
= − y (ϕ̂i,p+1/2 − ϕ̂i,p−1/2 ) ,
hi,j
(4.26)
onde foram denidas as grandezas
58
y
τ̄i,p
ademais, temos que
∫
1
= y
hi,j
yp+1/2
τ̂i (y)dy,
para
τi = ϕi
ou
yp−1/2 ≤ y ≤ yp+1/2
incógnitas na grade grossa exterior
m=M
temos
e
interna na
para
p=1
Γy
Λ.
xi−1/2 ≤ x ≤ xi+1/2 .
Especicamente para
x
˜x
J˜M
+1/2,j ≡ Ji+1/2,j ,
de fuga transversal na direção
(4.27)
p−1/2
Destacamos que algumas incógnitas da grade interna na
para
Jiy ,
x
x, F̃i,j
.
Γx
m=1
são coincidentes com
temos
y
y
Jˆi,1/2
≡ Jˆi,j−1/2
Analogamente, na outra direção, incógnitas da grade
Λ.
Sendo especíco,
e para
p=P
temos
y
ˆy
Jˆi,P
+1/2 ≡ Ji,j+1/2 ,
estas magnitudes
estão relacionadas com os termos de fuga transversal na direção
No problema dependente de
y
y , F̂i,j
.
x para cada nodo da grade interna Γx , as eq.(4.22) e (4.23)
representam um sistema de duas equações lineares e algébricas em quatro incógnitas:
x
J¯m,j
, ϕ̃m+1/2,j
e
grade interna
Γy ,
ϕ̄xm,j .
e
estas magnitudes estão relacionadas com os termos
são coincidentes com incógnitas na grade grossa exterior
temos
x
x
J˜1/2,j
≡ J˜i−1/2,j
De forma semelhante, no problema dependente de
y
x
J˜m+1/2,j
,
para cada nodo da
as eq.(4.25) e (4.26) representam um sistema de duas equações lineares e
algébricas em quatro incógnitas:
y
y
y
Jˆi,p+1/2
, J¯i,p , ϕ̂i,p+1/2 e ϕ̄i,p .
A m de obter uma única solução
para estes problemas, precisamos incluir duas equações auxiliares em cada nodo, considerar
as condições de continuidade nas interfaces das regiões e introduzir as condições de contorno
nas fronteiras do domínio.
As equações auxiliares propostas no método CSG-SGF-CN são
semelhantes às utilizadas no método SND-CN e aparecem na forma
ϕ̄xm,j =
ϕ̄yi,p =
γi,j
(γij − 1) y
F̂i,j ,
(ϕ̃m+1/2,j + ϕ̃m−1/2,j ) +
2
ΣTij
(4.28)
γi,j ˜x
x
x
J¯m,j
=
(Jm+1/2,j + J˜m−1/2,j
),
2
(4.29)
λi,j
(λij − 1) x
(ϕ̂i,p+1/2 + ϕ̂i,p−1/2 ) +
F̃i,j ,
2
ΣTij
(4.30)
λi,j ˆy
y
y
J¯i,p
=
(J
+ Jˆi,p−1/2
).
2 i,p+1/2
As eq.
(4.28) e (4.29) correspondem ao problema dependente de
correspondem ao problema dependente de
y.
(4.31)
x;
e as eq.
(4.28) e 4.29
Os parâmetros das equações auxiliares
γi,j
e
λi,j
devem ser determinados de forma que preservem no interior de cada nodo as soluções gerais
analíticas das equações de difusão integradas transversalmente dentro de cada célula externa
da grade externa.
Estas soluções estão dadas pelas Eqs.
(4.16 e 4.17) para a direção
x,
e
59
pelas eq. (4.18 e 4.19) para a direção
y.
Detalhes sobre o cálculo dos parâmetros
γi,j
e
λi,j
podem ser encontrados em [Barros et al. 2003].
As equações (4.22, 4.23, 4.28 e 4.29) em conjunto com as equações de continuidade
e as condições de contorno do lado esquerdo e direito do domínio, formam um sistema de
equações lineares e algébricas possível e determinado para o problema dependente de
x.
Ana-
logamente, as equações (4.25, 4.26, 4.30 e 4.31) em conjunto com as equações de continuidade
e as condições de contorno nas arestas superior e inferior do domínio, formam um sistema
de equações lineares e algébricas possível e determinado para o problema dependente de
y.
Destacamos que estes problemas podem ser resolvidos separadamente, entretanto eles estão
acoplados pelos termos de fuga transversal. Na próxima seção obtemos um esquema numérico
para implementar computacionalmente o método CSG-SGF-CN.
4.2 Esquema Numérico para o Método CSG-SND-CN
As equações (4.22), (4.23), (4.28) e (4.29) formam um sistema linear para o problema de
difusão unidimensional na direção
x.
As equações (4.25), (4.26), (4.30) e (4.31) formam um
sistema linear para o problema unidimensional na direção
são acoplados pelos termos de fuga transversal
y
(F̂i,j
ou
y.
x
F̃i,j
)
Estes problemas unidimensionais
através das arestas da célula
Λi,j
da malha grossa. Para convergir uma solução numérica do problema de difusão em geometria
X −Y,
usamos um esquema iterativo explícito de direções alternadas.
Para chegarmos no esquema de direções alternadas devemos simplicar os problemas
unidimensionais de cada direção. Iniciamos esse processo eliminando as correntes neutrônicas
presentes nos sistemas.
Para eliminarmos as correntes
x
J˜m±1/2,j
no problema da direção
x,
substituímos as equações auxiliares (4.28, 4.29) nas equações de balanço (4.22, 4.23) e os
resultados aparecem como
x
x
)=−
− J˜m−1/2,j
(J˜m+1/2,j
hxi,j ΣTij γi,j
y
(ϕ̃m+1/2,j + ϕ̃m−1/2,j ) − γi,j hxi,j F̂i,j
,
2
2D
x
x
)=−
+ J˜m−1/2,j
(J˜m+1/2,j
(ϕ̃m+1/2,j − ϕ̃m−1/2,j ) .
γi,j hxi,j
(4.32)
(4.33)
Em seguida somamos as equações (4.32) e (4.33) obtendo a seguinte expressão para a corrente
na aresta direita do nodo
60
]
)[ h Σ γ
Dij
xm Ti,j i,j
+
+
= ϕ̃m−1/2,j −
4
hxm γi,j
)
] (
(
)[ h Σ γ
Dij
hxm γij
xm Ti,j i,j
ϕ̃m+1/2,j −
−
γi,j −
F̂yi,j .
4
hx m
2
x
J˜m+1/2,j
(
(4.34)
Subtraindo as equações (4.32) e (4.33) obtemos a seguinte expressão para a corrente na aresta
esquerda do nodo
]
Dij
= ϕ̃m−1/2,j
+
+
4
hxm γi,j
] (
)
(
) [h Σ γ
Dij
hxm γij
xm Ti,j i,j
ϕ̃m+1/2,j
−
γi,j +
F̂yi,j .
4
hxm
2
x
J˜m−1/2,j
(
) [h
xm ΣTi,j γi,j
(4.35)
Neste ponto devemos aplicar as condições de continuidade da corrente nas interfaces
entre os nodos da grade interna
Γx .
Para isso consideramos
m = m+1
na eq.
(4.35) e
igualamos com a eq. (4.34). O resultado aparece na forma
[
]
hxij Σaij γij
Dij
ϕ̃m−1/2,j
−
+
4
hxij γij
]
[
hxij Σaij γij
2Dij
ϕ̃m+1/2,j
−
+
2
hxij γij
]
[
hxij Σaij γij
Dij
ϕ̃m+3/2,j
−
=
4
hxij γij
{
[
]
[
]
hxij νΣfij γij
hxij νΣfij γij
1
ϕ̃m−1/2
+ ϕ̃m+1/2,j
+
kef
4
2
[
]}
hxij νΣfij γij
− hxij γij F̂yij .
ϕ̃m+3/2,j
4
A eq.
(4.36) pode ser utilizada em cada aresta de nodo
grade grossa
Λij .
m ± 1/2
(4.36)
no interior da célula da
Para as arestas correspondentes às condições de contorno à esquerda e à
direita, devemos considerar as equações (3.11) e (3.12) respectivamente. Obtemos a equação
correspondente à aresta do contorno esquerdo considerando
m=1
e
i=1
na eq. (4.35) e
igualando o resultado à condição de contorno (3.11). O resultado aparece como
61
[
]
hx1j Σa1j γ1j
D1j
ϕ̃1/2,j αL +
+
+
4
hx1j γ1j
[
]
hx1j Σa1j γ1j
D1j
ϕ̃3/2,j
−
=
4
hx1j γ1j
{
[
]
hx1j νΣf1j γij
1
+
ϕ̃1/2,j
kef
4
]}
[
hx γ1j
hxij νΣfij γ1j
− 1j F̂y1j .
ϕ̃3/2,j
4
2
(4.37)
Obtemos a equação correspondente à aresta do contorno direito considerando
e
i = I
na eq.
m=M
(4.34) e igualando o resultado à condição de contorno (3.12), o resultado
aparece como
[
]
hxIj ΣaIj γIj
DIj
ϕ̃M −1/2,j
−
+
4
hxIj γIj
]
[
hxIj ΣaIj γIj
DIj
+
=
ϕ̃M +1/2,j αR +
4
hxIj γIj
{
[
]
hxIj νΣfIj γIj
1
ϕ̃M −1/2,j
+
kef
4
[
]}
hxIj νΣfIj γIj
hx γIj
ϕ̃M +1/2,j
− Ij F̂y1j
4
2
Os problemas unidimensionais na direção
malha grossa
Λ,
para
j = 1 : J.
contorno esquerdo (4.37),
M −1
x
(4.38)
deverão ser resolvidos em cada linha da
Em cada linha teremos uma equação para a aresta do
equações para as arestas interiores (4.36), e uma equação
para a aresta do contorno esquerdo (4.38). Os sistemas para cada linha
j
podem ser escritos
em notação matricial na forma
⃗x =
Ax ϕ
onde
Ax
1
Qx + F y ,
kef
são matrizes quadradas tridiagonais,
nas arestas verticais dos nodos (ϕ̃m±1/2,j ),
Qx
⃗x
ϕ
(4.39)
é o vetor de incógnitas contendo os uxos
contém os termos de fonte de ssão e
vetor contendo os termos de fuga transversal na direção
y
Fy
é um
em nível de região. As incógnitas
das camadas horizontais correspondentes aos problemas na direção
x
são mostradas na g.
4.3a.
Para eliminarmos as correntes neutrônicas
Jˆi,p±1/2 dos problemas na direção y seguimos
um procedimento semelhante. Primeiro substituímos as equações auxiliares (4.30, 4.31) nas
62
Figura 4.3: Incógnitas dos sistemas (a) Problema unidimensional dependente de
unidimensional dependente de
x,
(b) Problema
y.
equações de balanço (4.25, 4.26), tendo como resultado
y
y
(Jˆi,p+1/2
− Jˆi,p−1/2
)=−
hyi,j ΣTij λi,j
x
(ϕ̂i,p+1/2 + ϕ̂i,p−1/2 ) − λi,j hyi,j F̃i,j
,
2
2D
(Jˆi,p+1/2 + Jˆi,p−1/2 ) = −
(ϕ̂i,p+1/2 − ϕ̂i,p−1/2 ) .
λi,j hyi,j
(4.40)
(4.41)
Segundo, somamos as equações (4.40) e (4.41) obtendo a seguinte expressão para a corrente
na aresta superior do nodo
)
hyij ΣTij λij (
y
ˆ
ϕ̂i,p+1/2 + ϕ̂i,p−1/2 −
Ji,p+1/2 = −
4
) h λij
Dij (
y
ϕ̂i,p+1/2 − ϕ̂i,p−1/2 − ij
F̃xij .
hyij λij
2
(4.42)
Terceiro, subtraímos as equações (4.40) e (4.41) e obtemos a seguinte expressão para a corrente
na aresta inferior do nodo
)
hyij ΣTij λij (
y
ˆ
ϕ̂i,p+1/2 + ϕ̂i,p−1/2 −
Ji,p−1/2 =
4
) h λij
Dij (
y
ϕ̂i,p+1/2 − ϕ̂i,p−1/2 + ij
F̃xij .
hyij λij
2
(4.43)
63
Quarto, aplicamos as condições de continuidade da corrente nas interfaces entre os
nodos da grade interna
Γy .
Para isso consideramos
p=p+1
na eq. (4.43) e igualamos com
a eq. (4.42), o resultado aparece na forma
[
]
hyij Σaij λij
Dij
ϕ̂i,p−1/2
−
+
4
hyij λij
[
]
hyij Σaij λij
2Dij
+
ϕ̂i,p+1/2
+
2
hyij λij
[
]
hyij Σaij λij
Dij
−
ϕ̂i,p+3/2
=
4
hyij λij
[
{
]
hyij νΣfij λij
1
ϕ̂i,p−1/2
+
kef
4
[
]
hyij νΣfij λij
ϕ̂i,p+1/2
+
2
]}
[
hyij νΣfij λij
− hyij λij F̃xij .
ϕ̂i,p+3/2
4
A eq. (4.44) pode ser utilizada em cada aresta vertical de nodo
da grade grossa
(4.44)
p ± 1/2
no interior da célula
Λij .
Quinto e último passo, obter as equações para às arestas correspondentes as condições
de contorno superior e inferior. Para isso devemos considerar as equações (3.13) e (3.14) respectivamente. Obtemos a equação correspondente à aresta do contorno inferior considerando
p = 1 e j = 1 na eq.
(4.43) e igualando o resultado à condição de contorno (3.13). O resultado
aparece como
[
]
[
]
hyi1 Σai1 λi1
hyi1 Σai1 λi1
Di1
Di1
+ ϕ̂i,3/2
=
ϕ̂i,1/2 αB +
+
−
4
hyi1 λi1
4
hyi1 λi1
{
[
]
[
]} (
)
1
hyi1 νΣfi1 λi1
hyi1 νΣfi1 λi1
hyi1 λi1
ϕ̂i,1/2
+ ϕ̂i,3/2
−
F̃xi1 .
kef
4
4
2
Obtemos a equação correspondente à aresta do contorno superior considerando
e
j =J
(4.45)
p=P
na eq. (4.42) e igualando o resultado à condição de contorno (3.14). O resultado
aparece como
[
]
[
]
hyiJ ΣaiJ λiJ
DiJ
hyiJ ΣaiJ λiJ
DiJ
ϕ̂i,P −1/2
+
+ ϕ̂i,P +1/2 αT
−
=
4
hyiJ λiJ
4
hyiJ λiJ
{
[
]
[
]} (
)
hyiJ νΣfiJ λiJ
1
hyiJ νΣfiJ λiJ
hyiJ λiJ
ϕ̂i,P −1/2
+ ϕ̂i,P +1/2
−
F̃xiJ
kef
4
4
2
(4.46)
64
Os problemas unidimensionais na direção
da malha grossa
Λ,
para
contorno inferior (4.45),
i = 1 : I.
P −1
y
deverão ser resolvidos em cada coluna
Em cada coluna teremos uma equação para a aresta do
equações para as arestas horizontais interiores (4.44), e uma
equação para a aresta do contorno superior (4.46). Os sistemas para cada coluna
i podem ser
escritos em notação matricial na forma
⃗y =
Ay ϕ
onde
Ay
ϕ̂i,p±1/2
1
Qy + F x ,
kef
são matrizes quadradas tridiagonais, o vetor
que são as incógnitas do sistema,
os termos de fuga transversal na direção
Qy
x
⃗y
ϕ
(4.47)
contém os uxos escalares médios
contém os termos da fonte de ssão, e
Fx
contém
sobre células externas. As incógnitas das camadas
verticais dos sistemas são mostradas na gura 4.3b.
Os sistemas de equações (4.39) e (4.47) podem ser utilizados para calcularmos os uxos
nas arestas dos nodos para as direções
x
e
y
respectivamente. Estes sistemas estão acoplados
pelos termos de fuga transversal em nível de região dados pelas equações (4.13) e (4.15). Os
termos de fuga transversal são calculados em função das correntes neutrônicas nas arestas das
células externas
Λij .
Entretanto, as correntes neutrônicas foram eliminadas de nossos sistemas
de equações, precisamos então, expressar as correntes neutrônicas em função dos uxos.
Para obtermos uma nova expressão para as fugas na direção
m = 1 na eq.
(4.35) e
transversal na direção
m=M
x
x, F̃ijx ,
na eq. (4.34), e substituímos os resultados na denição de fuga
(4.15), obtemos então
[
]
[
]
hxij λij F̂yij + Hxij F̃xij = ϕ̃M +1/2,j + ϕ̃1/2,j χxij + ϕ̃M −1/2,j + ϕ̃3/2,j θijx
Para obtermos uma nova expressão para as fugas na direção
p=1
na eq. (4.43) e
transversal na direção
p=P
y
consideramos
y , F̂ijy ,
(4.48)
consideramos
na eq. (4.42), e substituímos os resultados na denição de fuga
(4.13), obtemos então
[
]
[
]
y
Hyij F̂yij + hyij F̃xij = ϕ̂i,P +1/2 + ϕ̂i,1/2 χyij + ϕ̂i,P −1/2 + ϕ̂i,3/2 θij
(4.49)
Nas eq. (4.48) e (4.49) consideramos as denições
χuij
huij ΣTij ∆uij
Dij
−
=−
4
huij ∆uij
(4.50)
65
u
θij
=−
onde
u=x
ou
y,
e
∆uij = γij
ou
huij ΣTij ∆uij
Dij
+
4
huij ∆uij
(4.51)
λij .
As equações (4.48) e (4.49) formam um sistema de duas equações lineares com duas
incógnitas:
F̃ijx
e
F̂ijy .
x
=
F̂i,j
Resolvemos o sistema de equações e obtemos
]
1 { y[ x
x
(ϕ̃m−1/2,j + ϕ̃3/2,j ) −
Hj χi,j (ϕ̃m+1/2,j + ϕ̃1/2,j ) + θi,j
ωij
[
]}
x
hxij χyi,j (ϕ̂i,P +1/2 + ϕ̂i,1/2 ) + θi,j
(ϕ̂i,P −1/2 + ϕ̂i,3/2 ) ,
(4.52)
]
1 { x[ y
y
(ϕ̂i,p−1/2 + ϕ̂i,3/2 ) −
Hi χi,j (ϕ̂i,p+1/2 + ϕ̂i,1/2 ) + θi,j
ωij
[
]}
y
y
x
hij χi,j (ϕ̃M +1/2,j + ϕ̃1/2,j ) + θi,j (ϕ̃M −1/2,j + ϕ̃3/2,j ) ,
(4.53)
y
F̂i,j
=
onde denimos:
ωi,j =
Hijx Hijx
1
.
− hxij γij hyij λij
(4.54)
Um esquema iterativo para convergir à solução numérica do problema bidimensional
é composto de iterações externas e internas. Nas iterações externas, novas estimativas para
o fator de multiplicação efetivo e para os termos de fonte de ssão são obtidas. Nas iterações
internas, usamos um esquema iterativo explícito de direções alternadas, segundo os seguintes
passos:
1. Primeiro usamos uma estimativa inicial dos uxos para calcularmos a fuga transversal
na direção
y
y , (F̂i,j
)
Eq. (4.53), na grade externa grossa
Λ.
2. Em seguida, calculamos o uxo escalar médio nas arestas verticais,
a eq. (4.39) sobre a grade espacial na
Γx
na direção
3. Então usamos o uxo calculado na grade na
as fuga na direção
x
)
x, (F̃i,j
Γx
(ϕ̃m±1/2,j ), resolvendo
x.
para obter uma nova estimativa para
eq. (4.52), sobre a grade externa grossa.
4. Com as fugas atualizadas calculamos o uxo médio nas arestas horizontais,
sobre a grade espacial interna na direção
y.
(ϕ̂i,p±1/2 ),
66
5. Os valores obtidos para
transversal na direção
y
(ϕ̂i,p±1/2 ),
são usados para gerar novas estimativas para a fuga
na grade espacial grossa e retornar ao passo 2.
As aproximações que consideramos no método CSG-SND-CN, onde os termos de fuga
transversal são aproximados ao longo das arestas das células de malha grossa, são mais rudimentares do que as aproximações do método clássico SND-CN, onde os termos de fuga
transversal são aproximados ao longo das arestas dos nodos. No entanto, esperamos que isso
não implique uma perda signicativa de precisão para determinar grandezas integrais não
localizadas, como o fator de multiplicação efetivo.
Resultados numéricos do método CSG-
SND-CN que ilustram esta hipóteses são oferecidos no próximo capítulo.
5 RESULTADOS NUMÉRICOS E DISCUSSÃO
Neste capítulo mostramos os resultados de diferentes experimentos numéricos que visam avaliar a precisão e eciência computacional do método CSG-SND-SGF. Os resultados
numéricos gerados pelo método CSG-SND-CN são comparados com aproximações analíticas e
resultados numéricos obtidos pelos métodos SND-CN e DD. Para este m, foram construídas
três rotinas em MATLAB que implementam os algoritmos de solução numérica de cada um
destes métodos.
Sendo o objetivo do método CSG-SND-SGF calcular grandezas integrais escolhemos
calcular o fator de multiplicação efetivo (kef ) em quatro problemas modelos típicos.
Dois
problemas com núcleo de reatores homogêneos e dois problemas com núcleos de reatores
heterogêneos.
Primeiro, mostramos os resultados para um reator homogêneo innito.
Em
seguida para um núcleo nito com as mesmas propriedades físicas. Segundo mostramos os
resultados um problema heterogêneo com condições de simetria e para um núcleo heterogêneo
de quatro regiões.
5.1 Reator Homogêneo Innito
O primeiro problema modelo resolvido foi um reator homogêneo innito, para modelar
este reator consideramos um domínio quadrado com
180 cm
reexivas nas quatro fronteiras do domínio, ou seja,
αL = αR = αT = αB = 0
3.7). Os dados nucleares considerados para o problema foram
e
de lado e condições de contorno
nas eq. (3.4-
D = 1.000 cm, Σa = 0.021 cm−1
νΣf = 0.022 cm−1 .
Resolvemos este problema modelo utilizando os três métodos numéricos descritos neste
trabalho e consideramos
ções externas) e
10−5
10−6
para o critério de convergência do autovalor dominante (itera-
para o critério de convergência da fuga de nêutrons nos esquemas de
direções alternadas (iterações internas). As rotinas numéricas foram executadas para diferentes tamanhos da grade espacial que variam de grades nas a grades grossas. Os resultados
numéricos para o fator de multiplicação efetivo e o numero de iterações externas necessárias
para a convergência são mostrados na tabela 5.1. Como solução de referência, consideramos
a seguinte aproximação analítica oferecida em [Duderstadt e Hamilton 1976]
νΣf
kef ∼
≈ 1, 047619 .
= k∞ =
Σa
(5.1)
68
Tabela 5.1: Fator de multiplicação efetivo e número de iterações externas do problema modelo de
reator innito.
Malha
90 × 90
a
1 × 90 × 90b
36 × 36
1 × 36 × 36
18 × 18
1 × 18 × 18
9 × 9
1 × 9 × 9
6 × 6
1 × 6 × 6
3 × 3
1 × 3 × 3
Método
kef
δk %c
Iterações
DD
1,047602
1,6273E-03
45
SND-CN
1,047620
9,0909E-05
40
CSG-SND-CN
1,047620
7,3727E-05
5
DD
1,047594
2,3909E-03
16
SND-CN
1,047613
5,7727E-04
31
CSG-SND-CN
1,047622
2,3409E-04
7
DD
1,047590
2,7727E-03
4
SND-CN
1,047611
7,6818E-04
14
CSG-SND-CN
1,047624
4,8227E-04
10
DD
1,0477587
3,0591E-03
3
SND-CN
1,047607
1,1500E-03
14
CSG-SND-CN
1,047628
8,2495E-04
17
DD
1,047514
1,0027E-02
4
SND-CN
1,047608
1,0545E-03
11
CSG-SND-CN
1,047632
1,2459E-03
20
DD
1,047486
1,2700E-02
14
SND-CN
1,047601
1,7227E-03
10
CSG-SND-CN
1,047635
1,4750E-03
31
a
grade espacial de I J, indica I células na direção x e J células na direção y
b
grades espacial composta de R I J, R células na grade externa, I nodos internos na direção
x, e J nodos internas na direção y
c
desvio percentual relativo com relação a solução de referencia
A partir da análise da tabela 5.1, podemos concluir que o método CSG-SND-CN gera
resultados precisos para o autovalor dominante no problema modelo de reator innito.
desvio relativo para todas as grades espaciais consideradas foi inferior a
0.1%.
O
O que pode
ser raticado se observamos a gura (5.1), que plota o desvio relativo percentual no cálculo
do
kef
nas grades consideradas.
69
0,015
DD
SND-CN
CSG-SND-CN
Desvio % relativo
0,010
0,005
0,000
90 x 90
36 x 36
18 x 18
9 x 9
6 x 6
3 x 3
Grade Espacial de I X J
Figura 5.1: Relação entre as grades espaciais e o desvio relativo para um reator innito.
Na gura (5.1) podemos perceber que nosso método apresenta desvios relativos percentuais bem próximos do método SND-CN. Os valores numéricos gerados pelo método CSGSND-CN são muito semelhantes aos valores numéricos gerados pelo método SND-CN. Neste
problema modelo, os termos de fuga transversal são zero, então as aproximações de ambos os
métodos são equivalentes. O que justica a forte correspondência entre os resultados obtidos.
Por enquanto, os resultados para o número de iterações externas não permitem colocar
nenhuma hipótese.
No método CSG-SND-CN o número de iterações diminuem quando as
grades se tornam mais nas. No entanto, no método SND-CN o número de iterações aumenta
quando se tornam mais nas as grades espaciais.
5.2 Reator Homogêneo Finito
O problema do reator homogêneo nito apresenta a mesma geometria e propriedades
físico-neutrônicas do domínio descrito na seção anterior. Consideramos condições de contorno
de uxo zero nas quatro fronteiras do domínio, ou seja,
αL = αR = αT = αB = ∞.
Resolvemos este problema modelo com os métodos DD, SND-CN e CSG-SND-CN
70
utilizando as rotinas computacionais desenvolvidas. Consideramos
vergência das iterações externas e
10−5
10−6
no critério de con-
no critério de convergência das iterações internas. Os
resultados numéricos obtidos para diferentes grades espaciais, que variam de uma grade na
de
90 × 90
com nodos quadrados de 2
nodos quadrados de 60
cm de lado,
cm
de lado até uma grade muito grossa de
3×3
com
são apresentados na tabela 5.2. Como solução de referên-
cia para o autovalor dominante consideramos a seguinte aproximação analítica, oferecida em
[Duderstadt e Hamilton 1976]
kef ∼
=[
(
Σa +
νΣf
)2
π
Lx
(
+
π
Ly
)2 ] ≈ 1, 018083 .
(5.2)
Como podemos observar na tabela 5.2, os valores gerados pelo método CSG-SND-CN
para o autovalor dominante neste segundo problema modelo tem uma boa precisão numérica;
em todas as grades espaciais consideradas os desvios relativos foram menores que
0.1%.
Um
gráco do desvio relativo percentual no cálculo do autovalor dominante dos três métodos é
oferecido na gura (5.2). Quando comparamos o CSG-SND-CN com o método espectro nodal
convencional SND-CN, a precisão numérica tem a mesma ordem em cálculos de malha grossa.
No entanto, a precisão numérica do método SND-CN é superior à precisão numérica do método
CSG-SND-CN para grades nas. Este comportamento é esperado, porque a aproximação das
fugas transversais em nível de região feita no método CSG-SND-CN , é mais grosseira que a
aproximação das fugas transversais em nível de nodo feita no método SND-CN.
71
Tabela 5.2: Fator de multiplicação efetivo e número de iterações externas para o problema modelo
do reator homogêneo nito.
Malha
90 × 90
a
1 × 90 × 90
36 × 36
b
1 × 36 × 36
18 × 18
1 × 18 × 18
9 × 9
1 × 9 × 9
6 × 6
1 × 6 × 6
3 × 3
1 × 3 × 3
Método
kef
δk %c
Iterações
DD
1,018086
2,7028E-04
74
SND-CN
1,018080
3,1907E-04
75
CSG-SND-CN
1,019033
9,3259E-02
91
DD
1,018088
4,6672E-04
75
SND-CN
1,018082
1,4226E-04
76
CSG-SND-CN
1,019032
9,3190E-02
90
DD
1,018147
6,2619E-03
77
SND-CN
1,018145
6,1107E-03
80
CSG-SND-CN
1,019033
9,3259E-02
88
DD
1,018364
2,7576E-02
80
SND-CN
1,018359
2,7114E-02
101
CSG-SND-CN
1,019032
9,3190E-02
80
DD
1,018722
6,2741E-02
84
SND-CN
1,018683
5,8910E-02
190
CSG-SND-CN
1,019031
9,3092E-02
77
DD
1,020619
2,4907E-01
63
SND-CN
1,018329
2,4139E-02
227
CSG-SND-CN
1,019032
9,3190E-02
76
72
0,5
DD
0,4
SND-CN
Desvio % relativo
CSG-SND-CN
0,3
0,2
0,1
0,0
90 x 90
36 x 36
18 x 18
9 x 9
6 x 6
3 x 3
Grade Espacial de I X J
Figura 5.2: Desvio relativo percentual no cálculo do autovalor dominante para o reator homogêneo
nito.
Na gura (5.3) plotamos o número de iterações externas dos métodos CSG-SND-CN
e SND-CN para as grades espaciais consideradas.
Observamos, que o número de iterações
do método CSG-SND-CN diminue suavemente quando as grades espaciais internas tornam-se
mais grossas. Por outro lado, o número de iterações do método SND-CN aumentam fortemente
quando aumentamos o tamanho dos nodos. Esse fato sugere que o método CSG-SND-CN é
mais eciente que o método SND-CN, em cálculos de malha grossa.
73
240
220
SND-CN
200
SCG-SND-CN
Iterações
180
160
140
120
100
80
60
90 x 90
36 x 36
18 x 18
9 x 9
6 x 6
3 x 3
Grade Espacial de I X J
Figura 5.3: Número de iterações externas segundo a grade espacial para o reator homogêneo nito.
Uma questão relevante é que o presente método espectronodal de grades compostas mostrou-se menos sensível a erros de discretização espacial do que o método SND-CN.
Outra característica do método CSG-SND-CN é que quando a grade espacial interna se
torna mais na, os resultados numéricos não se tornam mais precisos.
Esta característica
negativa foi observada para o método CSG-SGF-CN na formulação de ordenadas discretas
[Dominguez et al. 2009]. Este comportamento do método CSG-SND-CN é uma consequência da aproximação dos termos de fuga transversal nas arestas da malha externa (grossa).
Em outras palavras, a diminuição da largura das células espaciais
Γx
ou
Γy
tem uma fraca
inuência na precisão dos resultados.
5.3 Reator Heterogêneo Simétrico
O terceiro problema modelo é um núcleo de reator com condições de simetria, a geometria do problema e as condições de contorno são mostradas na gura (5.4). Observamos
condições de contorno reexivas nas fronteiras superior, inferior e esquerda; e condição de
contorno de uxo zero na fronteira direita. As propriedades físicas utilizadas são oferecidas
74
na tabela 5.3.
Figura 5.4: Geometria e condições de contorno para o problema de reator heterogêneo simétrico.
Tabela 5.3: Propriedades físicas para o problema de reator heterogêneo simétrico.
Material
Combustível I
Combustível II
D(cm) Σa (cm−1 ) νΣf (cm−1 )
1, 333333 0.033000
0, 035000
1, 333333 0.022000
0, 020000
Utilizando os três métodos descritos neste trabalho foram realizados diversos experimentos numéricos com o núcleo heterogêneo simétrico.
Os experimentos numéricos consi-
deraram diferentes grades espaciais, incluindo grades nas, médias e grossas.
No caso das
grades compostas utilizadas no método CSG-SND-CN na grade externa foi considerada uma
célula por região, ou seja, temos uma célula para a região de combustível I e uma célula para
a região de combustível II. Ademais foi considerado um critério de convergência de
10−6
no
desvio relativo do autovalor dominante (iterações externas), e a norma máxima entre as fugas
transversais de duas iterações sucessivas não pode superar
tados para o
kef
10−5
(iterações internas). Os resul-
e o número de iterações externas é mostrado na tabela 5.4. Para o cálculo do
desvio relativo consideramos como solução de referencia
métodos DD e SND-CN para uma grade na de
kef = 1, 043616
180 × 180
valor obtido com os
nodos. Para facilitar a analise dos
resultados, plotamos nas guras 5.5 e 5.6 o desvio relativo percentual e o numero de iterações
externas para cada uma das grades analisadas.
75
Tabela 5.4: Fator de multiplicação efetivo e numero de iterações externas para o problema heterogêneo
simétrico.
Malha
180 × 180
90 × 90
45 × 45
30 × 30
15 × 15
9 × 9
6 × 6
3 × 3
Método
kef
δk %c
Iterações
DD
1,043616
0,0000E+00
76
SND-CN
1,043616
0,0000E+00
69
CSG-SND-CN
1,043658
4,0245E-03
57
DD
1,043615
9,5821E-05
76
SND-CN
1,043616
0,0000E+00
69
CSG-SND-CN
1,043657
3,9286E-03
51
DD
1,043614
1,9164E-04
76
SND-CN
1,043618
1,9164E-04
72
CSG-SND-CN
1,043657
3,9286E-03
42
DD
1,043612
3,8328E-04
76
SND-CN
1,043620
3,8328E-04
77
CSG-SND-CN
1,043660
4,2161E-03
34
DD
1,043501
1,1019E-02
76
SND-CN
1,043627
1,0540E-03
95
CSG-SND-CN
1,043660
4,2161E-03
31
DD
1,043054
5,3851E-02
75
SND-CN
1,043662
4,4078E-03
110
CSG-SND-CN
1,043661
4,3119E-03
29
DD
1,042224
1,3338E-01
69
SND-CN
1,043702
8,2406E-03
131
CSG-SND-CN
1,043659
4,1203E-03
28
DD
1,040321
3,1573E-01
51
SND-CN
1,043946
3,1621E-02
142
CSG-SND-CN
1,043661
4,3119E-03
27
76
DD
0,035
SND-CN
CSG-SND-CN
0,030
Desvio % Relativo
0,025
0,020
0,015
0,010
0,005
0,000
-0,005
180 x 18090 x 90 45 x 45 30 x 30
15x15
9 x 9
6 x 6
3 x 3
Grade Espacial (I X J)
Figura 5.5: Desvio relativo percentual no cálculo do autovalor dominante para o problema modelo
do reator heterogêneo simétrico.
Da análise da tabela 5.4 e da gura 5.5 podemos constatar que o método CSG-SND-CN
obtém resultados numéricos satisfatórios para o autovalor dominante no problema heterogêneo
simétrico. Para todas as grades analisadas o desvio relativo obtido com o método CSG-SNDCN não supera o 0.1%.
Ademais enquanto a precisão dos resultados do método SND-CN
se deteriora para grades muito grossas, a precisão dos resultados no método CSG-SND-CN
praticamente não varia com as mudanças na grade.
Ao analisarmos a gura 5.6 observamos novamente que o número de iterações externas
do método CSG-SND-CN diminui suavemente para grades grossas. Enquanto o numero de
iterações externas no método SND-CN aumenta rapidamente em grades grossas. Este fato
sugere que o método CSG-SND-CN apresenta um melhor desempenho computacional que o
método SND-CN em cálculos de grades grossas obtendo resultados mais precisos com um
menor numero de iterações.
77
DD
140
SND-CN
SCG-SND-CN
Nº de Iterações
120
100
80
60
40
20
180 x 18090 x 90 45 x 45 30 x 30
15x15
9 x 9
6 x 6
3 x 3
Grade Espacial (I X J)
Figura 5.6: Numero de iterações externas para o problema modelo do reator heterogêneo simétrico.
5.4 Reator Heterogêneo de Quatro Regiões
O último problema modelo analisado é um reator heterogêneo de quatro regiões, resolvemos um quarto do núcleo do reator, a geometria do problema e as condições de contorno
são mostradas na gura 5.7. As propriedades físico-neutrônicas de cada um dos quatro tipos
de combustível aparecem na tabela 5.6.
Tabela 5.5: Propriedades físicas para o problema do reator de quatro regiões.
Material
Combustível I
Combustível II
Combustível III
Combustível IV
D(cm) Σa (cm−1 ) νΣf (cm−1 )
1, 200
0, 025
0, 026
1, 500
0.024
0, 025
1, 700
0.030
0, 032
1, 000
0, 021
0, 022
Utilizando as rotinas criadas para os métodos DD, SND-CN e CSG-SND-CN resolvemos o problema do reator de quatro regiões utilizando diferentes grades espaciais. Da mesma
forma que no problema modelo anterior, no método CSG-SND-CN a grade externa considera
uma célula espacial para cada região. Novamente consideramos um critério de convergência
78
Figura 5.7: Geometria e condições de contorno para o problema do reator heterogêneo de quatro
regiões.
de
10−6
nas iterações externas referente ao autovalor, e um critério de convergência de
nas iterações internas referente ao esquema de direções alternadas.
10−5
Como solução de refe-
rência consideramos os valores obtidos pelos métodos DD e SND-CN para uma grade na de
100×100, isto é kef = 1, 018459.
Os resultados numéricos obtidos para o autovalor dominante
e o número de iterações externas são mostrados na tabela 5.6.
A partir dos valores mostrados na tabela 5.6, construímos dois grácos, que mostram
a variação do desvio percentual relativo com relação à grossura da malha, que é a gura 5.8,
e o que mostra a variação do número de iterações em cada método com relação também a
grade espacial, gura 5.9.
79
Tabela 5.6: Fator de multiplicação efetivo e número de iterações externas para o problema heterogêneo
de quatro regiões.
Malha
100 × 100
50 × 50
20 × 20
10 × 10
8 × 8
4 × 4
Método
Kef
δK %c
Iteração
DD
1,018459
0,0000E+00
82
SND-CN
1,018459
0,0000E+00
92
CSG-SND-CN
1,018628
1,6594E-02
63
DD
1,018072
3,7999E-02
82
SND-CN
1,018451
7,8550E-04
96
CSG-SND-CN
1,018628
1,6594E-02
58
DD
1,013193
5,1706E-01
80
SND-CN
1,018526
6,5786E-03
99
CSG-SND-CN
1,018629
1,6692E-02
47
DD
1,003622
1,4568E+00
78
SND-CN
1,018746
2,8180E-02
114
CSG-SND-CN
1,018631
1,6888E-02
34
DD
1,001629
1,6525E+00
76
SND-CN
1,018906
4,3890E-02
132
CSG-SND-CN
1,018622
1,6005E-02
28
DD
0,974320
4,3339E+00
62
SND-CN
1,018853
3,8686E-02
151
CSG-SND-CN
1,018621
1,5906E-02
22
Da análise da tabela 5.6 e da gura 5.8 podemos armar que o método CSG-SNDCN calcula o autovalor dominante deste problema modelo com boa precisão numérica, para
todas as grades analisadas os desvios percentuais não superam o 0,2%. Ao compararmos o
método de grades compostas CSG-SND-CN com o método espectro nodal convencional SNDCN novamente observamos maior precisão do método de grades compostas em cálculos de
grades grossas. Também, neste problema modelo o método de grades compostas mostrou-se
menos sensível as mudanças nas dimensões dos nodos da grade espacial.
80
0,5
Desvio % Relativo
0,4
0,3
DD
SND-CN
CSG-SND-CN
0,2
0,1
0,0
-0,1
100 x 100
50 x 50
20 x 20
10 x 10
8 x 8
4 x 4
Grade Espacial (I X J)
Figura 5.8: Desvio relativo percentual no cálculo do autovalor dominante para o problema modelo
do reator heterogêneo de quatro regiões.
81
160
DD
SND-CN
140
CSG-SND-CN
Nº de Iterações
120
100
80
60
40
20
100 x 100
50 x 50
20 x 20
10 x 10
8 x 8
4 x 4
Grade Espacial (I X J)
Figura 5.9: Número de iterações externas para o problema modelo do reator heterogêneo de quatro
regiões.
Ao observarmos o gráco da gura 5.9 que mostra o número de iterações com relação
as dimensões da grade espacial, constatamos o mesmo comportamento que nos problemas
modelos anteriores. O método CSG-SND-CN apresenta uma diminuição no número de iterações enquanto engrossamos a grade espacial, ao contrário do método SND-CN que apresenta
um incremento do número de iterações.
6 CONCLUSÕES E TRABALHOS FUTUROS
Um dos principais desaos atuais da humanidade é a produção de energia com menor
impacto ao meio ambiente; especicamente minimizando a contribuição para o efeito estufa
e o aquecimento global. Mesmo após o acidente de Fukushima, a energia nuclear é uma das
alternativas promissoras para enfrentar este desao. Porém, dentre os entraves à expansão da
energia nuclear, está a melhoraria das ferramentas de modelagem e simulação computacional,
que garantam uma operação segura das usinas nucleares.
O presente trabalho propõe um novo método espectronodal de grades múltiplas, para
resolver problemas de autovalor de difusão em geometria cartesiana bidimensional e aproximação de uma velocidade. O método proposto utiliza uma aproximação constante para os
termos de fuga transversal a nível de região, e diminui o custo algébrico e computacional
dos método nodais convencionais, mantendo a precisão em cálculos de malhas grossas.
Os
fundamentos matemáticos do método foram apresentados, mostrando a discretização em dois
estágios: (i) malha externa grossa a nível de região e (ii) malha interna na a nível de nodos.
Deste modo, foram obtidos os sistemas de equações lineares e algébricos, que combinados,
formam um esquema explícito de direções alternadas para obter os uxos médios nos nodos,
e o fator efetivo de multiplicação.
O esquema numérico obtido foi implementado em rotinas computacionais do software
MATLAB, que permitem resolver problemas homogêneos e heterogêneos. As resoluções de
problemas modelos típicos permitem armar, que o método CSG-SND-CN é apropriado para
calcular grandezas integrais como o coeciente de multiplicação efetivo.
Nos experimentos
numéricos realizados, os erros para autovalor dominante, inclusive para grades muito grossas,
não ultrapassaram um ponto percentual.
Ao compararmos o CSG-SND-CN com o método espectronodal convencional (SNDCN), em relação à precisão dos resultados, constatamos:
•
Para as grades nas, o método SND-CN oferece resultados mais precisos que o método CSG-SND-CN. Este comportamento se justica, porque no método CSG-SND-CN
as fugas transversais são aproximadas por uma função constante em nível de região,
entretanto, no método SND-CN as fugas transversais são aproximadas por um valor
constante em nível de nodo.
Consequentemente o método CSG-SND-CN inclui uma
aproximação mais grosseira que o método SND-CN.
•
Para grades grossas e no cálculo de grandezas integrais, o método CSG-SND-CN ofe-
83
rece resultados mais precisos que o método nodal convencional. Este comportamento
se justica, porque o método de grades compostas mostrou-se menos sensível aos erros
de discretização espacial, ou seja, mudanças na grade espacial tinham pouca inuência na precisão dos resultados. Um elemento negativo do método CSG-SND-CN é que
um renamento nas malhas internas não reete uma melhoria na precisão dos cálculos, característica que foi observada por pesquisas anteriores com métodos de grades
compostas.
Ao compararmos o método proposto neste trabalho (CSG-SND-CN) com um método
espectronodal convencional (SND-CN), em relação ao custo computacional; podemos armar
que em cálculos de malha grossa, o método de grades compostas apresenta melhor desempenho
computacional. Esta armação é fundamentada no menor número de iterações externas apresentadas pelo método CSG-SND-CN, para a mesma grade espacial. Outro elemento relevante
é a diminuição do número de iterações que ocorre no método CSG-SND-CN ao diminuir o número de nodos na discretização. Ao contrário, se engrossarmos a malha no método SND-CN,
observamos um incremento no número de iterações.
O número de iterações não pode ser utilizado como critério absoluto sobre o desempenho computacional, entretanto, se compararmos as implementações dos métodos CSG-SNDCN e SND-CN, podemos observar que as iterações externas de ambos os métodos seguem
as mesmas etapas: (i) cálculo das fugas na direção
direção
x,
(iii) cálculo das fugas na direção
x,
y,
(ii) cálculo das grandezas nodais na
(iv) cálculo das grandezas nodais na direção
y,
(v) atualização da fonte de ssão e obtenção de uma nova estimativa para o autovalor, e (vi)
vericação do critério de convergência. Neste sentido, podemos armar que os custos computacionais destas iterações são equivalentes, nos permitindo utilizar o número de iterações
como critério de comparação. Em relação ao consumo de memória, as implementações dos
métodos também são equivalentes.
Outro elemento importante de comparação do método CSG-SND-CN com relação ao
método SND-CN é a diminuição do custo algébrico, para obter as equações discretizadas do
esquema numérico.
Os motivos desta diminuição estão no uso da grade espacial composta
e na discretização em dois estágios, que permitem dividir o problema bidimensional em dois
unidimensionais, que podem ser tratados separadamente. No caso do problema de autovalor
de difusão, a simplicação do trabalho algébrico não é tão signicativa quanto o que ocorre
em problemas de ordenadas discretas com fonte xa.
O futuro desdobramento desta pesquisa é a tradução das rotinas desenvolvidas em
MATLAB para a linguagem C ou FORTRAN, que permitem a avaliação do tempo de pro-
84
cessamento e geram códigos computacionais mais robustos e ecientes. Novos experimentos
numéricos que envolvam o cálculo de outras grandezas integrais como taxa de reação, potência
elétrica gerada ou fugas do núcleo devem ser realizados.
Em relação ao futuro dos métodos de grades compostas em problemas de autovalor de
difusão bidimensional, deve-se utilizar outras aproximações do tipo linear ou exponencial para
aproximar as fugas transversais em nível de região. Aproximações mais robustas devem conduzir a resultados numéricos mais precisos, sem aumentar signicativamente a complexidade
do método.
Referências Bibliográcas
[Azmy 1988]AZMY, Y. Y. Comparison of three approximations to the linear-linear nodal
transport method in weighted diamond-dierence form.
Nuclear Science and Engineering,
v. 100, p. 190200, 1988.
[Badruzzaman 1990]BADRUZZAMAN, A. Nodal methods in transport theory.
Advances in
Nuclear Science and Technology, v. 21, p. 293331, 1990.
[Barros e Filho 2007]BARROS, R. C.; FILHO, H. A. Progress in the formulation of the approximate albedo boundary conditions for one-speed x,y-geometry discrete ordinates and
diusion eigenvalue problems.
Progress in Nuclear Energy, v. 49, p. 161171, 2007.
[Barros et al. 2003]BARROS, R. C. et al. The application of spectral nodal methods to discrete ordinates and diusion problems in cartesian geometry for neutron multiplying systems.
Progress in Nuclear Energy, v. 42, p. 385426, 2003.
[Barros et al. 2004]BARROS, R. C. et al. Recent advances in spectral nodal methods for
numerically solving neutron diusion eigenvalue problems.
Transport Theory and Statistical
Physics, v. 33, p. 331346, 2004.
[Barros e Larsen 1992]BARROS, R. C.; LARSEN, E. W. A spectral nodal method for onegroup x,ygeometry discrete ordinates problems.
Nuclear Science and Engineering,
v. 111,
p. 3445, 1992.
[Barros e Silva 1999]BARROS, R. C.;
SILVA, F. C. Recent advances in spectral nodal
methods for x,ygeometry discrete ordinates deep penetration and eigenvalue problems.
Progress in Nuclear Energy, v. 35, p. 293331, 1999.
[Dominguez 2000]DOMINGUEZ, D. S.
Método Espectro Nodal de Difusión en dos Dimensi-
ones. Dissertação (Mestrado) Instituto Superior de Ciencias y Tecnologia Nuclear, 2000.
(
in Spanish).
[Dominguez e Barros 2007]DOMINGUEZ, D. S.; BARROS, R. C. The spectral green s function linearnodal method for one-speed x, y-geometry discrete ordinates deep penetration
problems.
Annals of Nuclear Energy, v. 34, p. 958966, 2007.
86
[Dominguez e Barros 2010]DOMINGUEZ, D. S.; BARROS, R. C. Spectral nodal method for
numerically solving two-energy group x,y geometry neutron diusion eigenvalue problems.
Int. J. Nuclear Energy Science and Technology, v. 5, n. 1, 2010.
[Dominguez et al. 2009]DOMINGUEZ, D. S. et al. Composite spatial grid spectral nodal
method for one-speed discrete ordinates deep penetration problems in x,y geometry.
Progress
in Nuclear Energy, 2009.
[Duderstadt e Hamilton 1976]DUDERSTADT, J.; HAMILTON, L.
Nuclear Reactor Analysis.
New York, USA: John Wiley and Sons, 1976.
[Goldemberg e Lucon 2007]GOLDEMBERG, J.; LUCON, O. Energia e meio ambiente no
brasil.
Estudos Avançados 21, v. 59, p. 720, 2007.
[Lamarsh e Baratta 2001]LAMARSH, J. R.; BARATTA, A. J.
Introduction to Nuclear Engi-
neering. Third. Upper Saddle Rider, New Jersey, USA: American Nuclear Society, 2001.
[Lawrence 1986]LAWRENCE, R. D. Progress in nodal methods for the solution of the neutron
diusion and transport equations.
Progress in Nuclear Energy, v. 17, p. 271301, 1986.
[Lewis e Miller 1993]LEWIS, E. E.; MILLER, W. F.
Computational Methods in Transport
Theory. rst. Illinois, USA: American Nuclear Society, 1993.
[Mendonça e Gutierez 2000]MENDONçA, M. J. C. de; GUTIEREZ, M. B. S.
O Efeito Estufa
e o Setor Energético Brasileiro. [S.l.], 2000.
[Stacey 2001]STACEY, W. M.
Nuclear Reactor Physics. NY, USA: John Wiley & Sons Inc,
2001.
[Zienkiewicz 1971]ZIENKIEWICZ, O. C.
2th. ed. NY, USA: McGraw Hill, 1971.
The Finite Element Methods in Engineering Science.