Anais do 12O Encontro de Iniciação Científica e Pós-Graduação do ITA – XII ENCITA / 2006
Instituto Tecnológico de Aeronáutica, São José dos Campos, SP, Brasil, Outubro, 16 a 19, 2006
Efeito das propriedades variáveis com o tempo em uma barra de um
reator nuclear
João Gilberto Furlan Rocha
Instituto Tecnológico de Aeronáutica - ITA/CTA
12228-900 – São José dos Campos, São Paulo, Brasil
Bolsista PIBIC-CNPq.
Endereço eletrônico: [email protected]
Edson Luis Zaparoli
Instituto Tecnológico de Aeronáutica - ITA/CTA
12228-900 – São José dos Campos, São Paulo, Brasil
Bolsista PIBIC-CNPq.
Endereço eletrônico: [email protected]
Resumo
Este projeto trata-se da simulação numérica da variação de temperatura em regimes permanente e
transiente no interior de uma barra de combustível de um reator nuclear. Esta é constituída de um cilindro
interior de UO2, revestido por uma camada de Zircaloy. As propriedades que foram utilizadas destes
elementos são: condutividade térmica, densidade, e calor específico; e cada uma delas é função da
temperatura. Dos dados obtidos, analisou-se a possibilidade de se implementar o mesmo problema
utilizando valores constantes para estas propriedades. As simulações foram feitas utilizando o software
Macsyma PDEase2D.
1. Introdução
A maior fonte conhecida de energia no universo material é, sem comparação, a energia nuclear, isto
é, a energia contida no núcleo de um átomo. Reatores nucleares são utilizados para realizar uma reação
nuclear de fissão com o objetivo de obter calor para alimentar uma termelétrica com ciclo a vapor. A
energia gerada pela fissão nuclear é transmitida ao fluido de resfriamento do núcleo do reator. A maioria
dos reatores nucleares utiliza barras cilíndricas agrupadas em feixes. As usinas nucleares são classificadas
em função do fluido que resfria o reator.
De acordo com [1], o PWR (reatores de água pressurizada) é o mais comum de todos os tipos de
reatores em uso no mundo, tanto para estações de energia comercial quanto para aplicações militares. O
combustível usado pelo PWR é o dióxido de urânio (UO2) em cantis de liga de zircônio. O UO2 é
compactado em pastilhas de 10 mm de diâmetro que têm a forma de um cilindro de cerca de 1 cm de
comprimento e de diâmetro, acomodadas no interior de varetas de 4,40 m de comprimento e 10,76 mm de
diâmetro de zircaloy (uma liga de zircônio), hermeticamente fechadas e que suportam temperaturas até
1.852 ºC (fig. 1). Cada vareta possui aproximadamente 2 kg de urânio. São necessárias 235 dessas varetas
para formar um conjunto, o elemento combustível.
Figura 1. Feixe de barras
1
Anais do 12O Encontro de Iniciação Científica e Pós-Graduação do ITA – XII ENCITA / 2006
Instituto Tecnológico de Aeronáutica, São José dos Campos, SP, Brasil, Outubro, 16 a 19, 2006
O objetivo deste projeto é analisar o comportamento da temperatura dentro destas barras, tanto em
regime permanente como em regime transitório, dada uma condição inicial, e também, analisar os efeitos
de se considerar o uso de propriedades dos materiais que as constituem constantes, comparando-os com
os resultados nos quais suas propriedades variam com a temperatura.
O CFD (Computational Fluid Dynamics) possibilita a análise computacional de problemas físicos
através de simulações [3]. Problemas que envolvem transferência de calor, fluido dinâmica, ou até mesmo
análise financeira podem ser implementados que utilizem o CFD. No caso deste projeto, utilizou-se o
software Macsyma PDEase2D. Independentemente do problema, gera-se uma malha que corresponde à
discretização do domínio do problema, e por métodos numéricos, calcula-se a solução em determinados
pontos, e efetua-se a interpolação dos resultados.
2. Formulação matemática e metodologia de solução
2.1. Definição matemática do problema
A equação que descreve a difusão de calor em três dimensões em um meio incompressível e sem
convecção é dada pela Eq. (1):
(1)
onde a solução é a temperatura T em função do tempo t e da posição, e:
k é a condutividade térmica do material (unidades no SI: W / m.K )
•
•
•
q& é a taxa de geração de calor no meio (unidades no SI: W / m3 )
c é o calor específico do material (unidades no SI: J / Kg °C )
p é a densidade do material (unidades no SI: Kg / m3 )
•
A taxa de geração de calor dentro do cilindro (fig. 2), utilizando-se de sua simetria cilíndrica, pode
⎡ ⎛ r ⎞2 ⎤
ser aproximada para, de acordo com [4]: g ( r ) = g 0 ⎢1 − ⎜ ⎟ ⎥ , onde b é o raio do elemento
⎣⎢ ⎝ b ⎠ ⎦⎥
combustível , r é a distância ao eixo central do cilindro, g 0 é constante.
Sejam: o cilindro combustível representado pelo índice 1; o Zircaloy, que reveste o cilindro interno,
pelo índice 2; a temperatura na face de contato entre os dois por TW; a temperatura do fluido exterior à
barra dada por T; e o coeficiente de transferência de calor entra a barra e o fluido dada por h.
a
K2
b
K1
g
2
Anais do 12O Encontro de Iniciação Científica e Pós-Graduação do ITA – XII ENCITA / 2006
Instituto Tecnológico de Aeronáutica, São José dos Campos, SP, Brasil, Outubro, 16 a 19, 2006
Figura 2. Elemento estudado.
As condições de contorno são dadas por:
d
(T1 ( 0 ) ) = 0
dr
T1 ( b ) = T2 ( b ) = TW
k2
d
(T2 ( a ) ) = h (T − T2 ( a ) )
dr
k1
d
d
T1 ( b ) ) = k2 (T2 ( a ) )
(
dr
dr
Observa-se que devido à simetria cilíndrica, e ao fato das barras serem consideradas longas, ocorre
variação da temperatura apenas na direção radial.
A condição inicial foi tomada como uma temperatura constante ao longo de toda a barra,
correspondente à temperatura do fluido de resfriamento.
2.2. Características físicas e propriedades dos materiais
Foi obtido de [1] que: o raio do cilindro interior vale 5mm; o raio da barra (combustível mais o
revestimento) vale 5,078mm. De [2] tem-se que o valor de g0 vale 4*107,8 W/m3 e que a temperatura do
fluido refrigerador vale 300ºC. O valor do coeficiente de transferência de calor h vale 450 W/m²ºC.
As propriedades dos materiais em cada região analisada são dadas por:
0< r < b
A partir de tabelas, fez-se o ajuste de curvas para temperaturas entre 273K e 1700K.
•
Densidade do UO2: p = -0,0003T + 11,062 (Mg/m3).A partir de valores da densidade de UO2
para determinadas temperaturas, realizou-se uma curva de ajuste destes valores, dada pela fig. 3.
densidade do UO2
11,4
densidade, Mg/m3
11,2
11
10,8
Seqüência1
10,6
Linear (Seqüência1)
10,4
y = -0,0003x + 11,062
R2 = 0,9973
10,2
10
0
500
1000
1500
2000
temperatura, K
Figura 3. Densidade do UO2 em função da temperatura.
3
Anais do 12O Encontro de Iniciação Científica e Pós-Graduação do ITA – XII ENCITA / 2006
Instituto Tecnológico de Aeronáutica, São José dos Campos, SP, Brasil, Outubro, 16 a 19, 2006
•
Condutividade térmica do UO2: K = 876,63T-0,7837 (W/m/K). Novamente, a partir de valores da
condutividade térmica do UO2, construiu-se uma curva de ajuste, dada pela fig. 4.
condutividade térmica, W/m/K
condutividade térmica do UO2
12
10
8
6
4
Seqüência1
Potência (Seqüência1)
2
0
0
500
1000
1500
2000
y = 876,63x -0,7837
R2 = 0,9985
temperatura, K
Figura 4. Condutividade térmica do UO2 em função da temperatura.
•
Calor específico [3]: 214,9 J/mol/K
b< r < a
Os valores dados para o Zircaloy foram obtidos de [2].
A condutividade térmica(W/m/K) do Zircaloy pode ser aproximada para temperaturas entre 300K e
1800K para:
• Condutividade térmica: K =12,767 – 5,4348x10-4 T + 8,9818x10-6 T 2
Para T < 1400 K, a densidade é dada por (Kg/m3):
• Densidade: p = 6595,2 – 0,1477T
Para 273 K < T < 1400 K, o calor específico é dado por (J/mol/K):
• Calor específico: c = 255,66 + 0,1024T
2.3. Metodologia de solução
De todos estas expressões apresentadas na secção anterior, e da equação que rege este problema,
percebe-se que o problema é descrito por uma equação diferencial parcial não-linear de segunda ordem,
com coeficiente variáveis, com uma variável (temperatura) dependente de duas variáveis independente:
raio (devido á simetria cilíndrica) e tempo.
Uma solução analítica para este caso é de difícil obtenção, assim, deve-se recorrer à solução
numérica. Existem três métodos tradicionais para a solução numérica de problemas de equações
diferenciais: Métodos das diferenças finitas, dos volumes finitos e dos elementos finitos. O software
utilizado para realizar as simulações foi o PDEase. Esse software utiliza o método de elementos finitos de
Galerkin com uma malha não-estruturada para a representação do domínio e discretização das equações
diferenciais parciais. Foi utilizado o refinamento automático de malha nas regiões de maior gradiente de
temperatura com o objetivo de reduzir o esforço computacional. Desta forma, o sistema original de
equações diferenciais parciais foi transformado em um sistema de equações algébricas, resolvido por
meio de algoritmos iterativos. A fig 5 mostra a malha que o software utiliza para determinar a solução.
Quanto maior a concentração de pontos, maior é o esforço computacional para determinar a solução.
Observa-se uma concentração maior de pontos na junção entre os dois materiais da barra.
4
Anais do 12O Encontro de Iniciação Científica e Pós-Graduação do ITA – XII ENCITA / 2006
Instituto Tecnológico de Aeronáutica, São José dos Campos, SP, Brasil, Outubro, 16 a 19, 2006
Figura 5. Malha utilizada na solução.
3. Análise de resultados
3.1. Regime permanente
Para o caso em que a condução ocorre em regime permanente, as curvas de nível da temperatura na
barra são dadas pela fig. 6. Observa-se que a temperatura no interior se mantém entre 982°C e 1381°C.
Um gráfico da distribuição de temperatura, tomando-se uma secção da barra, em função do raio, está dado
na fig. 7. Nota-se que há uma mudança brusca na temperatura entre a superfície da barra e o fluido
refrigerador.
Figura 6. Curvas de nível da temperatura em regime permanente.
5
Anais do 12O Encontro de Iniciação Científica e Pós-Graduação do ITA – XII ENCITA / 2006
Instituto Tecnológico de Aeronáutica, São José dos Campos, SP, Brasil, Outubro, 16 a 19, 2006
Figura 7. Temperatura em função da distância ao centro.
3.2. Regime transitório
Durante o regime transiente, até que se atinja o regime estacionário, a temperatura no elemento
combustível varia de 573K a 1654K. E no revestimento esta variação se dá de 573K a 1260K. A fig. 8
indica a temperatura em função do tempo para 4 pontos da barra: um deles no centro, outro com r = 0,2a,
outro com r = 0,5a e outro com r = 0,8a.
Estamos interessados em trabalhar com coeficientes constantes de densidade, capacidade térmica e
de condutividade. Para isso, tomaram-se os valores médios de cada um deles no intervalo em que ocorre a
variação de temperatura dada. Estes valores (no S.I.) estão dados na Tab. 1, e foram obtidos baseados nas
fórmulas da secção 2.2.
Tabela 1. valores médios das propriedades dos meios para o intervalo de variação de temperatura
Condutividade térmica
densidade
Capacidade térmica
UO2
4,33
10728
215
Zircaloy
20,87
12920
350
A fig. 9 indica a temperatura em função do tempo para 4 pontos da barra: um deles no centro, outro
com r = 0,2a, outro com r = 0,5a e outro com r = 0,8a. Neste caso, considerou-se a propriedades acima
constantes. Percebe-se que ocorrem mudanças no comportamento entre as duas soluções. Na solução da
figura 9, o aumento da temperatura se dá de uma maneira mais rápida que a outra solução, e atinge o
regime estacionário mais rapidamente, e a temperaturas maiores. Estas podem ser observadas também
através de curvas de níveis tomando-se, por exemplo, o instante t = 20s, no qual se observam diferenças
significativas entre as respostas. A fig. 11 mostra isso.
Figura 9: temperatura em função do tempo para quatro pontos distintos com coeficientes variáveis
Figura 10. temperatura em função do tempo para quatro pontos distintos com coeficientes constantes
6
Anais do 12O Encontro de Iniciação Científica e Pós-Graduação do ITA – XII ENCITA / 2006
Instituto Tecnológico de Aeronáutica, São José dos Campos, SP, Brasil, Outubro, 16 a 19, 2006
Figura 11. Observação das diferenças entre as respostas para propriedades variáveis e constantes, no
instante t = 20s.
4. Agradecimentos
Aos professores Edson Luis Zaparoli e Cláudia Regina de Andrade, que me orientaram neste projeto.
5. Referências
[1] http://www.energiatomica.hpg.ig.com.br/ (consultado em novembro e dezembro 2005)
[2] Hagrman, D. T., ed., “SCDAP/RELAP5/MOD 3.1 Code Manual: MATPRO - A Library of
Materials Properties for Light-Water-Reactor Accident Analysis, NUREG/CR-6150, EGG2720 Vol. 4 (1995).
[3] Meliska. C.F.D. Computational Fluid Dynamics
[4] Özisik, M. Necati. Heat transfer – A basic approach
[5] Perry, Robert H. Perry’s Chemical Engineers’ Handbook
[6] Shaw, R. K. Heat Transfer Equipment Design
7