Acessibilidade / Reportar erro

Determinação de Fontes de Nêutrons que Conduzem Sistemas Subcríticos a Distribuições Prescritas de Potência

RESUMO

Um sistema nuclear é dito subcrítico quando os eventos responsáveis pela remoção dos nêutrons (fugas pelos contornos estruturais do sistema e absorção) acontecem em maior intensidade que os eventos que promovem a produção destas partículas (fissão). Quando isto acontece o sistema não consegue manter um nível estável em relação à população de nêutrons e tende ao desligamento. Por outro lado, qualquer sistema subcrítico pode ser dirigido por uma ou mais fontes estacionárias de nêutrons. O objetivo deste trabalho é apresentar uma metodologia para determinar as intensidades das fontes internas estacionárias, isotrópicas e uniformes de nêutrons que devem ser inseridas em um sistema subcrítico para dirigi-lo a uma distribuição prescrita de potência. Para tanto, foi utilizada a equação unidimensional do transporte de nêutrons monoenergéticos com espalhamento isotrópico em meios multiplicativos, referenciada como equação física de transporte, e a equação que é adjunta a esta equação, ambas na formulação das ordenadas discretas (SN). O método de malha grossa matriz resposta adjunto (RM) foi aplicado para resolver numericamente as equações SN adjuntas. O método RM é um método da mesma classe do método matriz resposta utilizado na solução das equações SN físicas de transporte. Resultados numéricos para dois problemas-modelo típicos são apresentados para examinar a metodologia proposta.

Palavras-chave:
ordenadas discretas; problema adjunto; matriz resposta adjunto; sistemas subcríticos

ABSTRACT

A nuclear system is referred to as subcritical when neutron removal events (leakage through the boundaries of the system and absorption) occur at a rate which is greater than the events that promote the production of these particles (fission). When this occurs the system cannot maintain a stable level of the population of neutrons and tends to shutdown. On the other hand, any subcritical system may be driven by one or more stationary sources of neutrons. The purpose of this work is to present a methodology to determine the intensities of the stationary and uniform sources of neutrons that must be inserted in a subcritical system to drive it to generate a prescribed power distribution. To achieve this goal, we use the time-independent, slab-geometry, one-speed neutron transport equation as a mathematical model in multiplying media, referred to as forward problem, and the equation which is adjoint to it, both in discrete ordinates (SN) formulation. The coarsemesh adjoint response matrix (RM) method is applied to numerically solve the adjoint SN problems. The RM method is a companion method of the response matrix method used to solve forward SN equations. Numerical results for two typical model problems are presented to verify the offered methodology.

Keywords:
discrete ordinates; adjoint problem; adjoint response matrix; subcritical systems

1 INTRODUÇÃO

Os reatores nucleares apresentam-se hoje como uma alternativa altamente viável no que diz respeito a fontes de energia limpa. Fusão e fissão, reações nucleares que liberam uma significativa quantidade de energia, não produzem dióxido de carbono (CO2). Junto a isto, o processo de implantação das usinas termonucleares não altera em grande escala o ecossistema local onde a usina será situada. Estes fatores, dentre outros, contribuem para que a energia nuclear possua um interessante potencial sob a ótica tecnológica.

Dentre o desenvolvimento de novas tecnologias na área nuclear surge o reator ADS (Accelerator Driven Systems). O ADS é um dispositivo que possui um núcleo subcrítico, isto é, os eventos que promovem a remoção de nêutrons acontecem em maior quantidade do que os eventos que promovem a produção destas partículas. Diante disso, este reator necessita do auxílio de fonte(s) externa(s) e uniforme(s) de nêutrons para manter o seu funcionamento 1212 S.A. Pereira. “Um conceito alternativo de um reator híbrido (conjunto subcrítico acoplado com acelerador)”. Ph.D. thesis, IPEN, Instituto de Pesquisas Energéticas e Nucleares, São Carlos, SP (2002).. Esta característica gera um intrínseco ganho na perspectiva de segurança, em relação aos atuais reatores nucleares comerciais, devido à possibilidade de desligamento da fonte caso alguma irregularidade de funcionamento ocorra. Como o sistema é subcrítico, o desligamento da(s) fonte(s) externa(s) de nêutrons induzirá o desligamento quase imediato do reator 1111 H. Nifenecker, O. Meplan & S. David. “Accelerator Driven Subcritical Reactors”. Institute of Physics Publishing, U.K. (2003)..

Neste contexto, o presente trabalho tem como objetivo apresentar uma metodologia que permite determinar as intensidades das fontes internas estacionárias e uniformes de nêutrons que devem ser inseridas em um sistema subcrítico para dirigi-lo a uma distribuição prescrita de potência. Para tanto, considera-se o problema unidimensional do transporte de nêutrons monoenergéticos, com espalhamento isotrópico e condições de contorno reflexiva ou vácuo como problema físico de transporte a ser tratado. Além disso, faz-se referência ao problema adjunto como um problema de valor de contorno onde os fluxos angulares adjuntos emergentes do sistema são nulos e as fontes adjuntas definidas no interior do domínio são conhecidas. A clássica formulação das ordenadas discretas (SN) é aplicada nos problemas físicos e adjuntos, onde é utilizada a quadratura de Gauss-Legendre 99 E.E. Lewis & W.F. Miller. “Computacional methods of neutrons transport”. American Nuclear Soc, New York (1993). para aproximar os termos integrais existentes nas equações que modelam matematicamente estes problemas.

Para resolver numericamente as equações SN adjuntas foi aplicado um método de malha grossa, inserido na classe dos métodos espectronodais, denominado matriz resposta adjunto (RM). Fundamentado no método matriz resposta 1515 O.P. Silva. “Um Método de Matriz Resposta para Cálculos de Transporte Multigrupos de Energia na Formulação das Ordenadas Discretas em Meios Não-Multiplitcativos.”. Ph.D. thesis, Universidade do Estado do Rio de Janeiro, Nova Friburgo, RJ (2018)., utilizado na solução das equações SN físicas de transporte, o método RM† faz uso de um conjunto completo de soluções elementares das equações SN adjuntas em cada região constituinte do domínio para construir o esquema iterativo adjunto de inversão por região (RBI) que converge resultados numéricos para os fluxos angulares adjuntos nas interfaces das regiões constituintes do domínio, que são livres de erro de truncamento espacial.

A seguir um breve resumo do conteúdo deste trabalho é exibido. Na seção 2 é apresentado o modelo matemático adotado, bem como o desenvolvimento de uma metodologia para determinar as fontes estacionárias e uniformes de nêutrons que devem ser inseridas em um sistema subcrítico para dirigi-lo a uma distribuição prescrita de potência. Na seção 3 é descrito o método RM que será utilizado na obtenção das soluções dos problemas SN adjuntos. Na seção 4 os resultados numéricos de experimentos realizados em dois problemas-modelo são apresentados e na seção 5 uma breve discussão sobre o trabalho com sugestões para trabalhos futuros é realizada.

2 MODELO MATEMÁTICO E A METODOLOGIA ADJUNTA

Dividimos esta seção em duas partes. Na seção 2.1, são apresentados os modelos matemáticos utilizados e na seção 2.2, a metodologia seguida para estimar as fontes de nêutrons.

2.1 Modelo matemático

As equações que modelam matematicamente os problemas unidimensionais físicos e adjuntos, para um sistema subcrítico estabilizado por fontes estacionárias e uniformes de nêutrons possuem respectivamente as seguintes formas:

μ x ψ x , μ + σ t x ψ x , μ = σ s x + v σ f x 2 - 1 1 ψ x , μ ' d μ ' + Q x (2.1)

- μ x ψ x , μ + σ t x ψ x , μ = σ s x + v σ f x 2 - 1 1 ψ x , μ ' d μ ' + Q x , (2.2)

onde por definição 44 D.G. Cacuci. “Handbook of nuclear engineering”. Springer, New York (2010).:

µ = cos(θ) = direção de migração dos nêutrons, onde θ é o angulo polar medido em relação ao eixo x;

ψ(x,µ) = fluxo angular de nêutrons;

ψ (x,µ) = fluxo angular adjunto;

σt (x) = seção de choque macroscópica total;

σs (x) = seção de choque macroscópica de espalhamento isotrópico;

σf (x) = seção de choque macroscópica de fissão;

ν = número médio de nêutrons liberados em um evento de fissão;

Q(x) = fonte interna estacionária e uniforme de nêutrons;

Q (x) = fonte adjunta.

As equações (2.1) e (2.2) aparecem na notação de operadores como

H ψ = Q (2.3)

H ψ = Q . (2.4)

Aqui são omitidas as variáveis independentes (x, µ) no fluxo angular de nêutrons e no fluxo angular adjunto e a variável independente (x) na fonte isotrópica de nêutrons e na fonte adjunta, apenas por questão de simplificação de notação. Claramente, os operadores lineares H e H , que aparecem nas equações (2.3) e (2.4), são respectivamente definidos como

H μ x · + σ t x · - σ s x + v σ f x 2 - 1 1 · d μ ' , (2.5)

H - μ x · + σ t x · - σ s x + v σ f x 2 - 1 1 · d μ ' , (2.6)

onde (·) é um símbolo de posição. Apesar de as equações (2.1) e (2.2) diferirem apenas no sinal do termo de migração, elas possuem significados diferentes. Na equação (2.1), Q(x) é a intensidade da(s) fonte(s) internas de nêutrons e ψ(x,µ) pode ser entendido como a população de nêutrons está distribuída no espaço de fase. Na equação (2.2), Q (x) é a fonte adjunta, uma quantidade arbitrária que não possui uma relação física com Q(x). A grandeza ψ (x,µ) pode ser entendida no contexto deste trabalho como a importância, ou provável contribuição que os nêutrons em um ponto do espaço de fase tem para a constituição de uma certa quantidade integral mensurável. Esta quantidade integral pode assumir um determinado significado físico se escolhermos Q (x) de modo adequado 88 J. Lewins. “Importance the adjoint function”. Pergamont Press, U.K. (1965).. São exemplos desta quantidade integral: a taxa de absorção dos nêutrons em uma certa região do domínio e a leitura de um detector de nêutrons posicionado em um ponto do sistema 1616 G.E. Sjoden . Deterministic adjoint transport applications for He-3 neutron detector design. Annals ofNuclear Energy , 29 (2002), 1055-1071.), (1313 K. Royston, A. Haghighat, W. Walters, C. Yi & G.E. Sjoden. Determination of current at a detector window using a hybrid adjoint function methodology. Progress in Nuclear Science and Technology, 4 (2014), 528-532.) e (66 J.P. Curbelo, O.P. da Silva, C.R. García & R.C. Barros. Generalization of Spectral Green’s Function nodal method for slab-geometry fixed-source adjoint transport problems in SN formulation. In “International Nuclear Atlantic Conference-INAC 2017”. Associação Brasileira de Energia Nuclear (ABEN), Belo Horizonte, Brazil (2017)..

Para obter as soluções das equações (2.1) e (2.2) é preciso definir as condições de contorno a serem consideradas em cada tipo de problema. As condições de contorno utilizadas em problemas físicos de transporte aplicados em física de reatores nucleares são convencionalmente de dois tipos: vácuo ou reflexiva. No caso da condição de contorno tipo vácuo, não há a possibilidade de nenhuma partícula entrar pelos contornos estruturais do sistema, isto é, o fluxo angular de nêutrons em todas as direções que incidem nos contornos do sistema é igual a zero. No caso da condição de contorno reflexiva, os valores dos fluxos angulares de nêutrons que incidem no sistema são iguais aos valores dos fluxos angulares de nêutrons que emergem do sistema. Esta condição de contorno é utilizada em problemas que possuem simetria física e geométrica dos materiais, onde neste caso os fluxos angulares de nêutrons emergentes no ponto de simetria do domínio são refletivos de volta, gerando assim os fluxos angulares incidentes neste ponto. Desta forma, ao invés de resolver o problema percorrendo o domínio espacial inteiro, é adotada a condição de contorno reflexiva, que permite resolver o mesmo problema de forma reduzida partindo do ponto de simetria do domínio. Assim, o esforço computacional gasto na solução de um problema é reduzido.

No caso do problema adjunto as condições de contorno tradicionalmente utilizadas são a de fluxo angular adjunto nulo nas direções emergentes dos contornos e a reflexiva. Na condição de contorno tipo fluxo adjunto emergente nulo, a importância dos nêutrons que emergem dos contornos do sistema é nula, pois como não há possibilidade de esses nêutrons retornarem ao sistema, eles não irão contribuir com a quantidade integral. Isto significa dizer que o fluxo angular adjunto, nas direções emergentes dos contornos do sistema, é igual a zero. A condição de contorno reflexiva nos problemas adjuntos de transporte é similar à do problema físico. Aqui cabe ressaltar que, assim como a equação adjunta é vinculada a equação física de transporte, as condições de contornos utilizadas nos problemas adjuntos são vinculadas as condições de contorno utilizadas nos problemas físicos. A mesma situação ocorre quando a condição de contorno adotada para o problema físico é a reflexiva.

2.2 Metodologia

Nesta seção apresenta-se uma metodologia para estimar as fontes estacionárias e uniformes que devem ser inseridas em um sistema subcrítico para dirigi-lo gerando uma distribuição prescrita de potência. Para tanto, é descrita a seguir uma expressão que será usada no desenvolvimento da metodologia

ψ , H ψ = ψ , H ψ . (2.7)

A equação (2.7) é conhecida na literatura como relação de reciprocidade 44 D.G. Cacuci. “Handbook of nuclear engineering”. Springer, New York (2010).. Aqui 〈 (·), (·)〉 é o produto interno entre funções e significa a integração em todo o domínio das variáveis independentes. Isto é,

f 1 , f 2 = 0 L - 1 1 f 1 x , μ f 2 x , μ d μ d x ,

onde f 1 e f 2 são duas funções definidas no mesmo espaço de fase de ψ e ψ e L é o comprimento do domínio espacial considerado neste trabalho (viz. Figura 1).

Figura 1:
Domínio unidimensional com k = 1 : K.

Neste ponto, faz-se a observação que a equação (2.7) só é válida na forma apresentada devido às condições de contorno adotadas neste trabalho para a solução das Eqs. (2.1) e (2.2). Logo, multiplicando a equação (2.4) por ψ, a equação (2.3) por ψ , aplicando o produto interno nas duas equações resultantes e considerando a equação (2.7) obtém-se

ψ , Q = ψ , Q . (2.8)

A variável Q deve ser escolhida de modo a transformar a quantidade integral do lado esquerdo da equação (2.8) em um valor de potência. Assim, temos uma relação entre o fluxo angular adjunto, a distribuição de fontes e o valor de potência gerado. Com esta escolha, o fluxo angular adjunto pode ser entendido como a importância, ou provável contribuição que os nêutrons em um ponto do espaço de fase terão para a geração desta potência. Definindo Q como

Q x = ε σ f x δ k x , (2.9)

onde ε representa a energia média liberada em uma reação de fissão, i.e., ε=200MeV e δkx=1,xRk 0,caso contrário, com R k representando uma dada região combustível do domínio, tem-se

ψ , ε σ f x δ k x = ψ k , Q . (2.10)

O índice usado no fluxo angular adjunto remete à escolha da fonte adjunta; portanto, ψk é o fluxo angular adjunto obtido através da equação (2.2) com fonte adjunta dada pela equação (2.9). Neste ponto, analisaremos os termos que constituem a equação (2.10). No termo à esquerda tem-se

ψ , ε σ f x δ k x = 0 L ε σ f x δ k x - 1 1 ψ x , µ ' d μ ' d x , (2.11)

onde a integral -11ψ(x,µ')dµ' é definida como fluxo escalar de nêutrons ϕ(x). A equação (2.11) pode, portanto, ser reescrita como

ψ , ε σ f x δ k x = 0 L ε σ f x δ k x ϕ x d x = P k , (2.12)

onde P k significa a potência gerada pela região R k . Em continuidade, no termo à direita da equação (2.10), tem-se

ψ k , Q = - 1 1 Q x - 1 1 ψ k x , µ ' d μ ' d x , (2.13)

onde a integral -11 ψk(x,μ')dμ' é definida como fluxo escalar adjunto ϕk(x). Portanto,

ψ k , Q = 0 L Q ( x ) ϕ k ( x ) d x . (2.14)

Igualando as equações (2.12) e (2.14) obtém-se

P k = 0 L Q ( x ) ϕ k ( x ) d x . (2.15)

Visando à simulação computacional deste problema discretiza-se a variável angular nas equações (2.1) e (2.2). O método de discretização utilizado é a clássica formulação das ordenadas discretas (SN), que consiste na discretização da variável angular µ em N direções e na utilização de conjuntos de quadraturas, que neste trabalho serão as quadraturas de Gauss-Legendre, para aproximar os termos integrais existentes nas equações que modelam os problemas físicos e adjuntos. Assim, em uma região homogeneizada R k do domínio apresentado na Figura 1, as equações (2.1) e (2.2) na formulação SN aparecem respectivamente como

μ m d d x ψ m ( x ) + σ t k ψ m ( x ) = σ s k + ν σ f k 2 n = 1 N ψ n ( x ) ω n + Q k , m = 1 : N , x k - 1 / 2 < x < x k + 1 / 2 , (2.16)

e

- μ m d d x ψ m ( x ) + σ t k ψ m ( x ) = σ s k + ν σ f k 2 n = 1 N ψ n ( x ) ω n + Q k , m = 1 : N , x k - 1 / 2 < x < x k + 1 / 2 . (2.17)

Aqui N é a ordem da quadratura, µ m as direções discretas de migração dos nêutrons e ω m são os pesos da quadratura. Por simplicidade é utilizado nas equações (2.16) e (2.17) as notações ψm(x)=ψ(x,μm) e ψm(x)=ψ(x,μm). Neste trabalho os parâmetros materiais, fontes de nêutrons e fontes adjuntas são considerados uniformes em relação à variável espacial no interior das regiões que compõem o domínio. Portanto, para xk-1/2<x<xk+1/2:σt(x)=σtk, σs(x)=σsk, σf(x)=σfk, Q(x)=Qk, Q(x)=Qk. Com as aproximações apresentadas acima, a equação (2.15) torna-se

P k = l = 1 K Q l x l - 1 / 2 x l + 1 / 2 ϕ k ( x ) d x . (2.18)

Multiplicando a equação (2.18) por hlhl, onde h l é a espessura da região R l e aplicando o conceito de valor médio de uma função, obtém-se

P k = l = 1 K Q l ϕ ¯ l , k h l , (2.19)

onde ϕ¯l,k é definido como fluxo escalar médio adjunto na região R l gerado pela fonte adjunta na região R k . Repetindo este processo e variando k=1:K, constrói-se o seguinte sistema matricial

P = L Q , (2.20)

onde P é um vetor coluna de ordem K, composto pelas potências prescritas; Q é um vetor coluna de ordem K, composto pelas fontes de nêutrons; e L é uma matriz quadrada de ordem K, definida como

L = h 1 ϕ ¯ 1 , 1 h 2 ϕ ¯ 2 , 1 h K ϕ ¯ K , 1 h 1 ϕ ¯ 1 , 2 h 2 ϕ ¯ 2 , 2 h K ϕ ¯ K , 2 h 1 ϕ ¯ 1 , K h 2 ϕ ¯ 2 , K h K ϕ ¯ K , K . (2.21)

Repetindo este processo para k=1:K faz-se necessário resolver K problemas adjuntos variando a fonte adjunta de acordo com a região k. Assim, as fontes de nêutrons que devem ser inseridas em um sistema subcrítico para estabilizá-lo e, com isso, gerar uma distribuição prescrita de potência são

Q = L - 1 P . (2.22)

3 O MÉTODO MATRIZ RESPOSTA ADJUNTO

Esta seção é dividida em duas partes. Na seção 3.1 são construídas, em cada região homogeneizada do sistema, um conjunto de soluções elementares para o problema adjunto e na seção 3.2 utilizam-se estas soluções para construir um conjunto específico de matrizes que permitem estimar os valores do fluxo angular adjunto nas interfaces de cada região homogeneizada do sistema.

3.1 Soluções gerais locais

A solução geral da equação (2.17) em uma região R k tem o seguinte formato 55 K.M. Case & P.F. Zweifel. “Linear transport theory”. Addison-Wesley Publishing, Massachusetts (1967).

ψ m , k ( x ) = ψ m h ( x ) + ψ k p , m = 1 : N , x k - 1 / 2 < x < x k + 1 / 2 , (3.1)

onde os superíndices h e p indicam respectivamente as componentes homogênea e particular da solução geral local. A ausência do subíndice m na componente particular da solução geral apresentada na equação (3.1) deve-se à característica isotrópica da fonte de nêutrons que consideramos neste trabalho.

Para determinar a componente homogênea da solução geral local, deve-se resolver a seguinte equação

- μ m d d x ψ m h ( x ) + σ t k ψ m h ( x ) = σ s k + ν σ f k 2 n = 1 N ψ n h ( x ) ω n , m = 1 : N , x k - 1 / 2 < x < x k + 1 / 2 . (3.2)

Seguindo os trabalhos 22 L.B. Barichello & C.E. Siewert. A Discrete Ordinates Solution For A Non-Grey Model With Complete Frequency Redistribution. Journal of Quantitative Spectroscopy and Radiative Transfer , 62 (1999), 665-675.),(66 J.P. Curbelo, O.P. da Silva, C.R. García & R.C. Barros. Generalization of Spectral Green’s Function nodal method for slab-geometry fixed-source adjoint transport problems in SN formulation. In “International Nuclear Atlantic Conference-INAC 2017”. Associação Brasileira de Energia Nuclear (ABEN), Belo Horizonte, Brazil (2017)., a componente homogênea da equação (3.2) possui a seguinte forma

ψ m , ϑ h ( x ) = a m ( ϑ ) e - ( x - γ ) ϑ , m = 1 : N , x k - 1 / 2 < x < x k + 1 / 2 , (3.3)

onde γ=xk-1/2,ϑ>0xk+1/2,ϑ<0.Na Eq. (3.3) é usada a estratégia do deslocamento na função exponencial como proposto na literatura 11 L.B. Barichello, R.D.M. Garcia & C.E. Siewert. A Spherical HarmonicsSolution For RadiativeTransfer Problems With Reecting Boundaries AndInternal Sources. Journal of Quantitative Spectroscopy and Radiative Transfer, 60(2) (1998), 247-260.),(22 L.B. Barichello & C.E. Siewert. A Discrete Ordinates Solution For A Non-Grey Model With Complete Frequency Redistribution. Journal of Quantitative Spectroscopy and Radiative Transfer , 62 (1999), 665-675.),(33 M. Benassi, R.M. Cotta & C.E. Siewert. The PN Method for RadiativeTransfer Problems with Reective Boundary Conditions. Journal of Quantitative Spectroscopy and Radiative Transfer , 30 (1983), 547-553.),(77 J.P. Curbelo , O.P. da Silva , C.R. García & R.C. Barros . “Integral Methods in Science and Engineering”, volume 2. Springer International Publishing, 1 ed. (2017), chapter Shifting Strategy in the Spectral Analysis for the Spectral Green s Function Nodal Method for Slab-Geometry Adjoint Transport Problems in the Discrete Ordinates Formulation, pp. 201-210.),(1717 R. Vargas & M.T. Vilhena. A closed-form solution for one-dimensional radiative condutive problem by the decomposition and LTSn methods. Journal of Quantitative Spectroscopy and Radiative Transfer , 61 (1999), 303-308.. Esta estratégia limita as N funções exponenciais da solução apresentada na Eq. (3.3) no intervalo (0, 1), evitando assim overflow exceptions na solução deste problema devido à aritmética finita computacional.

Substituindo a equação (3.3) na equação (3.2) obtém-se

n = 1 N σ s k + ν σ f k 2 μ m ω n - σ t k μ m δ m , n a n = 1 ϑ a m , (3.4)

onde δm,n=1,m=n0,caso contrário.Na sequência, variando m=1:N na equação (3.4), obtém-se um sistema matricial na forma de um problema de autovalor

X a = 1 ϑ a , (3.5)

onde a é um autovetor de N componentes correspondente a cada autovalor 1ϑ e X é uma matriz quadrada de ordem N. Resolvendo o problema apresentado na equação (3.5), obtêm-se N autovalores reais e distintos e seus respectivos autovetores. Os autovalores obtidos possuem a seguinte característica: são compostos por N/2 valores positivos e N/2 valores negativos. Neste trabalho os autovalores foram ordenados em uma forma crescente. A componente homogênea pode ser escrita como uma combinação linear das soluções propostas na equação (3.3), com os autovalores e seus respectivos autovetores dados pela solução do problema apresentado na equação (3.5). Logo

ψ m h ( x ) = = 1 N α ψ m , ϑ h ( x ) = = 1 N α a m ( ϑ ) e - ( x - γ ) ϑ , m = 1 : N , x k - 1 / 2 < x < x k + 1 / 2 , (3.6)

onde α 𝓁 são constantes arbitrárias.

O próximo passo é determinar a componente particular da solução geral local. Devido ao fato de a componente particular ser constante na região R k e independente da direção de propagação dos nêutrons, a sua determinação é uma tarefa mais simples. Substituindo ψkp na equação (2.17) tem-se

ψ k p = Q k σ t k - ( σ s k + ν σ f k ) . (3.7)

De posse das componentes homogênea e particular, a solução geral local aparece como

ψ m , k ( x ) = = 1 N α a m ( ϑ ) e - ( x - γ ) ϑ + Q k σ t k - ( σ s k + ν σ f k ) , m = 1 : N , x k - 1 / 2 < x < x k + 1 / 2 . (3.8)

Neste ponto, faz-se a ressalva, que este procedimento deve ser realizado em todas as regiões constituintes do sistema.

3.2 A matriz resposta adjunta

Tendo encontrado a solução geral local em todas as regiões do sistema, são determinados em cada região os coeficientes α ?? para o cálculo dos valores do fluxo angular adjunto nessas regiões. O método RM faz uso dos fluxos angulares adjuntos emergentes nas interfaces de cada região como forma de estimar os coeficientes α 𝓁 e criar matrizes que serão utilizadas para calcular os fluxos angulares adjuntos incidentes iterativamente, varrendo da esquerda para a direita μm<0 e da direita para esquerda μm>0 até que um critério de parada prescrito seja atingido.

Este procedimento é iniciado escrevendo a equação (3.8) na seguinte forma matricial

Ψ k ( x ) = A d i a g e - ( x - γ ) ϑ α + Ψ k p , x k - 1 / 2 < x < x k + 1 / 2 , (3.9)

onde Ψk(x) é um vetor coluna composto pelos fluxos angulares adjuntos em cada direção discreta m=1 : N; A é uma matriz quadrada de ordem N, composta pelos N autovetores obtidos através da solução da equação (3.5); diage-(x-γ)ϑ é uma matriz diagonal de ordem N, composta pelos elementos não-nulos e-(x-γ)ϑ, com =1 : N; α é um vetor coluna de ordem N, composto pelos coeficientes de expansão; e Ψkp é um vetor coluna descrito pela equação (3.7). A equação (3.9) pode ser reescrita como

Ψ 1 , k ( x ) Ψ 2 , k ( x ) Ψ N , k ( x ) = a 1 ( ϑ 1 ) e - ( x - γ ) ϑ 1 a 1 ( ϑ N ) e - ( x - γ ) ϑ N a 2 ( ϑ 1 ) e - ( x - γ ) ϑ 1 a 2 ( ϑ N ) e - ( x - γ ) ϑ N a N ( ϑ 1 ) e - ( x - γ ) ϑ 1 a N ( ϑ N ) e - ( x - γ ) ϑ N α 1 α 2 α N + Ψ k p . (3.10)

Para os fluxos angulares adjuntos presentes na equação (3.10), considerando μm<0 e μm>0, faz-se, respectivamente x=xk-1/2 e x=xk+1/2. Assim, a equação 3.10 assume a forma

Ψ μ m > 0 , k ( x k + 1 / 2 ) Ψ μ m < 0 , k ( x k - 1 / 2 ) = A 1 , 1 A 1 , 2 d i a g e - h k ϑ = ( N / 2 + 1 ) : N A 2 , 1 d i a g e - h k ϑ = 1 : N / 2 A 2 , 2 α + Ψ k p , (3.11)

onde Ψμm>0,k é um vetor coluna de ordem N/2, composto pelos fluxos angulares adjuntos com μm>0; Ψμm<0,k é um vetor coluna de ordem N/2, composto pelos fluxos angulares adjuntos com μm<0; e A1,1, A1,2, A2,1 e A2,2 são submatrizes quadradas de ordem N/2 que compõem a matriz A . Definindo

M = A 1 , 1 A 1 , 2 d i a g e - h k ϑ = ( N / 2 + 1 ) : N A 2 , 1 d i a g e - h k ϑ = 1 : N / 2 A 2 , 2 , (3.12)

pode-se expressar α a partir das equações (3.11) e (3.12). Logo,

α = M - 1 Ψ μ m > 0 , k ( x k + 1 / 2 ) Ψ μ m < 0 , k ( x k - 1 / 2 ) - Ψ k p . (3.13)

Na sequência, substituindo a equação (3.13) na (3.10) tem-se

Ψ μ m > 0 , k ( x ) Ψ μ m < 0 , k ( x ) = T ( x ) Ψ μ m > 0 , k ( x k + 1 / 2 ) Ψ μ m < 0 , k ( x k - 1 / 2 ) + [ I N - T ( x ) ] Ψ k p , (3.14)

onde I N é uma matriz identidade de ordem N e T(x) é uma matriz quadrada de ordem N denominada matriz resposta adjunta

T ( x ) = A d i a g e - ( x - γ ) ϑ M - 1 . (3.15)

Fazendo x=xk+1/2 na equação (3.14), encontra-se

Ψ μ m > 0 , k ( x k + 1 / 2 ) Ψ μ m < 0 , k ( x k + 1 / 2 ) = I N / 2 O N / 2 G - G + Ψ μ m > 0 , k ( x k + 1 / 2 ) Ψ μ m < 0 , k ( x k - 1 / 2 ) + O N / 2 O N / 2 - G - I N / 2 - G + Ψ k p . (3.16)

Na equação (3.16) O N/2 é a matriz nula de ordem N/2 e I N/2 é a matriz identidade de ordem N/2. As matrizes G+ e G- são matrizes quadradas de ordem N/2 que serão usadas para estimar os fluxos angulares adjuntos nas direções incidentes na interface à direita da região R k , isto é, os fluxos angulares adjuntos com μm<0 no ponto x k+1/2 . Então,

Ψ μ m < 0 , k ( x k + 1 / 2 ) = G + Ψ μ m < 0 , k ( x k - 1 / 2 ) + G - Ψ μ m > 0 , k ( x k + 1 / 2 ) + F k , (3.17)

onde Fk=(IN/2-G+)Ψinf,kp-G- Ψsup,kp. Ψsup,kp é a metade superior e Ψinf,kp é a metade inferior do vetor coluna Ψkp. Aqui vale a pena ressaltar que devido a aplicação da quadratura de GaussLegendre na aproximação dos termos integrais existentes nas equações SN físicas e adjuntas, temos necessariamente N como um número par. O mesmo procedimento deve ser realizado para estimar os fluxos angulares adjuntos nas direções incidentes na interface à esquerda da região R k , isto é, os fluxos angulares adjuntos com μm>0 no ponto x k−1/2 . Fazendo x=xk-1/2 na equação (3.14), obtemos

Ψ μ m > 0 , k ( x k - 1 / 2 ) = G + Ψ μ m > 0 , k ( x k + 1 / 2 ) + G - Ψ μ m < 0 , k ( x k - 1 / 2 ) + F k . (3.18)

Impondo condições de continuidade nas interfaces das regiões é feita uma varredura no domínio espacial da esquerda para a direita, através da equação (3.17) iniciando com as condições de contorno Ψμm<0(0) e depois é executada essa varredura da direita para a esquerda através da equação (3.18) iniciando com as condições de contorno Ψμm>0(L). Este esquema de varredura é realizado sucessivamente até que um critério de parada seja satisfeito. O critério de parada adotado neste trabalho requer que o desvio relativo entre duas estimativas consecutivas dos fluxos escalares adjuntos sejam menores ou iguais a ξ=1×10-5. Este esquema iterativo é denominado esquema iterativo adjunto de inversão por região (RBI ).

Um esquema iterativo similar é aplicado para resolver o problema físico e encontra-se descrito na referência 1414 D.M. Silva, E.J. Lydia, M.R. Guida, J.H. Zani & H.A. Filho. Analytical methods for computacional modeling of fixed-source slab-geometry discrete ordinates transport problems: response matrix and hybrid SN . Progress inNuclear Energy , 69 (2013), 77-84.. Uma vantagem na aplicação do método RM neste tipo de problema é que a fonte adjunta não faz parte do processo construtivo das matrizes G+ e G-. Portanto, se conservado os parâmetros materiais do problema, pode-se realizar várias simulações computacionais alterando a fonte adjunta, sem a necessidade de calcular as matrizes G +e G em cada simulação.

4 RESULTADOS NUMÉRICOS

Nesta seção são apresentados dois problemas-modelo unidimensionais e heterogêneos: o primeiro problema é composto por cinco regiões multiplicativas; e o segundo problema é constituído por três regiões não multiplicativas e duas multiplicativas.

No primeiro problema-modelo serão realizados quatro ensaios numéricos. Este problema é ilustrado na Figura 2, onde C.C.V indica que as condições de contorno utilizadas foram do tipo vácuo para o problema físico e do tipo fluxo adjunto nulo nas direções emergentes dos contornos para o problema adjunto. Os parâmetros materiais das três zonas que constituem esse domínio estão listados na Tabela 1. A quadratura escolhida foi a S4 de Gauss-Legendre e a potência total ser gerada neste problema-modelo foi considerada como 800 MW. O fator de multiplicação efetivo deste sistema é keff=0, 89157 que o insere no contexto dos reatores ADS 1111 H. Nifenecker, O. Meplan & S. David. “Accelerator Driven Subcritical Reactors”. Institute of Physics Publishing, U.K. (2003)..

Figura 2:
Problema-modelo 1.

Tabela 1:
Parâmetros materiais do problema-modelo 1.

No primeiro ensaio numérico é definida uma distribuição simétrica de potência, onde as regiões 1 e 5 geram 10% da potência total cada uma, as regiões 2 e 4 geram 20% da potência total cada uma e a região 3 gera 40% da potência total. Com essa distribuição de potência as fontes determinadas foram: Q1=Q5=1,49814 × 1016nêutrons/cm 3 s, Q2=Q4=1,86195 × 1016nêutrons/cm 3 s e Q3=1,31933 × 1017nêutrons/cm 3 s. Para confirmar se esta distribuição única de fontes guia o sistema subcrítico a distribuição de potência informada, foi simulado e resolvido um problema S4 físico com a distribuição de fontes determinada e calculada as potências geradas por cada região do domínio. Os valores das potências geradas em cada região neste caso foram: 80 MW para cada uma das regiões 1 e 5, 160 MW para cada uma das regiões 2 e 4, e 320 MW para a região 3.

No segundo ensaio numérico a distribuição de potência é definida do seguinte modo: as regiões 1 e 5 geram cada uma 5% da potência total; as regiões 2 e 4 geram cada uma 10% da potência total e a região 3 gera sozinha 70% da potência total. Com essa distribuição de potência as fontes geradas foram: Q1=Q5=9,31416 × 1015nêutrons/cm 3 s, Q2=Q4=-5,67535 × 1015nêutrons/cm 3 s e Q3=2,40531 × 1017nêutrons/cm 3 s. Do ponto de vista prático essa distribuição de potência é fisicamente absurda, pois não existe fonte negativa de nêutrons. Isto é um indicativo que a distribuição de potência deve ser mais suave para fazer sentido, pois para que se gere 70% da potência total na região 3, é necessário que uma grande quantidade de nêutrons esteja confinada nesta região. Esses nêutrons, quando não induzirem fissões poderão migrar para as regiões vizinhas que geram uma potência muito baixa (10%); portanto, a fonte negativa é um indicativo de que nêutrons devam ser retirados destas regiões para que se possa conseguir os 10% de potência, que é um valor muito baixo em comparação à região central. Para ilustrar esta assertiva, suavizamos a distribuição de potência no terceiro ensaio, onde agora é definida da seguinte forma: as regiões 1 e 5 geram cada uma 5% da potência total; as regiões 2 e 4 geram cada uma 15% da potência total e a região 3 gera 60% da potência total. Neste caso não ocorreu a determinação de fontes negativas, que foram de valores: Q1=Q5=7,94424 × 1015nêutrons/cm 3 s, Q2=Q4=5,49804 × 1015nêutrons/cm 3 s e Q3=2,03557 × 1017nêutrons/cm 3 s. Realizando o mesmo procedimento do primeiro ensaio para confirmar a validade numérica das fontes encontradas, tem-se que as potências geradas por região são: 40 MW para as regiões 1 e 5, 120 MW para as regiões 2 e 4, e 480 MW para a região 3.

Concluímos o primeiro problema-modelo com o quarto ensaio numérico, onde considera-se uma distribuição de potência uniforme; neste caso 20% da potência total gerada em cada região. As fontes geradas foram: Q1=Q5=3,17956 × 1016nêutrons/cm 3 s, Q2=Q4=2,25157 × 1016nêutrons/cm 3 s e Q3=6,26305 × 1016nêutrons/cm 3 s. Aqui, foi observado que distribuições uniformes de potência não implicam necessariamente em distribuições uniformes de fontes para dirigirem o sistema subcrítico. Esta característica se deve ao fato do domínio ser composto por regiões com espessuras e seções de choques macroscópicas diferentes.

No segundo problema-modelo foi realizado apenas um ensaio numérico. Este problema encontrase ilustrado na Figura 3 e os parâmetros materiais das três zonas que constituem este domínio estão listados na Tabela 2. Neste experimento foi utilizada a quadratura S4 de Gauss-Legendre e a potência total a ser gerada neste problema-modelo foi considerada, assim como no primeiro problema-modelo, como 800 MW. Com esta configuração o fator de multiplicação efetivo keff=0,80238, resultado que nos insere no contexto dos reatores ADS 1111 H. Nifenecker, O. Meplan & S. David. “Accelerator Driven Subcritical Reactors”. Institute of Physics Publishing, U.K. (2003)..

Figura 3:
Problema-modelo 2.

Tabela 2:
Parâmetros materiais do problema-modelo 2.

Neste ensaio numérico foi definida uma distribuição de potência, onde a região 2 gera 1% da potência total e a região 4 gera sozinha os 99% restantes. Com essa distribuição de potência as fontes uniformes geradas foram: Q1=Q3=Q5=0nêutrons/cm 3 s, Q2=3,24577 × 1015nêutrons/cm 3 s e Q4=3,21375 × 1017nêutrons/cm 3 s. Para confirmar se esta distribuição de fontes guia o sistema subcrítico a distribuição de potência informada foi solucionado um problema S4 físico utilizando as fontes determinadas e calculamos as potências geradas por cada região combustível do domínio. Os valores das potências geradas em cada região foram: 0 MW para as regiões 1, 3 e 5; 8 MW para a região 2 e 792 MW para a região 4. Pode ser observado que, apesar da distribuição de potência não ser nada suave, não foram gerados valores negativos para as fontes internas. Isso se deve a alta seção de choque de espalhamento da região 3, o que faz com que grande parcela dos nêutrons que migram para esta região sejam refletidos de volta para a região de origem. Sendo assim, a parcela de nêutrons que migram da região 4 e que chegam na região 2 é muito pequena. Então poucas fissões serão induzidas por esses nêutrons, o que gera um valor positivo de fonte interna necessária para esta região.

5 CONCLUSÕES

Foi apresentada e analisada neste trabalho uma metodologia para estimar as intensidades das fontes isotrópicas e uniformes de nêutrons que devem ser inseridas em um sistema subcrítico para dirigi-lo a uma distribuição prescrita de potência. Para tanto, foram apresentados na seção 2 os modelos matemáticos utilizados para os problemas físico e adjunto e descrito uma metodologia que permite estimar as fontes a partir da entrada das potências geradas por cada região combustível do domínio. Devido à complexidade da obtenção da solução analítica das equações que modelam os problemas físicos e adjuntos, foi utilizada a formulação SN para discretizar a variável angular µ em N direções e assim transformar as equações analíticas, que modelam os problemas deste trabalho, em sistemas de equações mais simples de serem resolvidos numericamente. Para resolver numericamente estes conjuntos de equações, foi detalhado na seção 3 o método matriz resposta adjunto que faz uso de um completo conjunto de soluções elementares em cada região homogeneizada do domínio. Com as condições de contorno e de continuidade, o método RM em conjunto com o esquema iterativo RBI estima os valores dos fluxos angulares adjuntos nas interfaces de cada região. Na seção 4 foram realizados 5 ensaios numéricos com o intuito de analisar a aplicabilidade da metodologia desenvolvida.

Para trabalhos futuros, sugerimos a inclusão da dependência energética, que leva em consideração as transferências de energia cinética entre nêutrons e núcleos-alvo nos eventos de espalhamento e fissão nuclear 1010 R.L. Murray. “Nuclear Energy”. Elsevier, New York (2009).. Outro ponto interessante a ser trabalhado, reside na extensão do modelo matemático adotado considerando uma geometria multidimensional cartesiana. Por fim, a realização de uma análise focada nos aspectos computacionais da metodologia é necessária e deve aguardar trabalhos futuros. Considerando principalmente problemas multidimensionais, o processo de construção da matriz-importância deve se tornar computacionalmente custoso, o que fomenta a necessidade da utilização de algoritmos paralelizados na sua construção.

AGRADECIMENTOS

Este trabalho foi financiado em parte pela Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES). Os autores gostariam de expressar seu reconhecimento ao Laboratório Multi-escala e Transporte de Partículas (LABTRAN) no IPRJ/UERJ devido ao suporte técnico prestado.

REFERÊNCIAS

  • 1
    L.B. Barichello, R.D.M. Garcia & C.E. Siewert. A Spherical HarmonicsSolution For RadiativeTransfer Problems With Reecting Boundaries AndInternal Sources. Journal of Quantitative Spectroscopy and Radiative Transfer, 60(2) (1998), 247-260.
  • 2
    L.B. Barichello & C.E. Siewert. A Discrete Ordinates Solution For A Non-Grey Model With Complete Frequency Redistribution. Journal of Quantitative Spectroscopy and Radiative Transfer , 62 (1999), 665-675.
  • 3
    M. Benassi, R.M. Cotta & C.E. Siewert. The PN Method for RadiativeTransfer Problems with Reective Boundary Conditions. Journal of Quantitative Spectroscopy and Radiative Transfer , 30 (1983), 547-553.
  • 4
    D.G. Cacuci. “Handbook of nuclear engineering”. Springer, New York (2010).
  • 5
    K.M. Case & P.F. Zweifel. “Linear transport theory”. Addison-Wesley Publishing, Massachusetts (1967).
  • 6
    J.P. Curbelo, O.P. da Silva, C.R. García & R.C. Barros. Generalization of Spectral Green’s Function nodal method for slab-geometry fixed-source adjoint transport problems in SN formulation. In “International Nuclear Atlantic Conference-INAC 2017”. Associação Brasileira de Energia Nuclear (ABEN), Belo Horizonte, Brazil (2017).
  • 7
    J.P. Curbelo , O.P. da Silva , C.R. García & R.C. Barros . “Integral Methods in Science and Engineering”, volume 2. Springer International Publishing, 1 ed. (2017), chapter Shifting Strategy in the Spectral Analysis for the Spectral Green s Function Nodal Method for Slab-Geometry Adjoint Transport Problems in the Discrete Ordinates Formulation, pp. 201-210.
  • 8
    J. Lewins. “Importance the adjoint function”. Pergamont Press, U.K. (1965).
  • 9
    E.E. Lewis & W.F. Miller. “Computacional methods of neutrons transport”. American Nuclear Soc, New York (1993).
  • 10
    R.L. Murray. “Nuclear Energy”. Elsevier, New York (2009).
  • 11
    H. Nifenecker, O. Meplan & S. David. “Accelerator Driven Subcritical Reactors”. Institute of Physics Publishing, U.K. (2003).
  • 12
    S.A. Pereira. “Um conceito alternativo de um reator híbrido (conjunto subcrítico acoplado com acelerador)”. Ph.D. thesis, IPEN, Instituto de Pesquisas Energéticas e Nucleares, São Carlos, SP (2002).
  • 13
    K. Royston, A. Haghighat, W. Walters, C. Yi & G.E. Sjoden. Determination of current at a detector window using a hybrid adjoint function methodology. Progress in Nuclear Science and Technology, 4 (2014), 528-532.
  • 14
    D.M. Silva, E.J. Lydia, M.R. Guida, J.H. Zani & H.A. Filho. Analytical methods for computacional modeling of fixed-source slab-geometry discrete ordinates transport problems: response matrix and hybrid SN . Progress inNuclear Energy , 69 (2013), 77-84.
  • 15
    O.P. Silva. “Um Método de Matriz Resposta para Cálculos de Transporte Multigrupos de Energia na Formulação das Ordenadas Discretas em Meios Não-Multiplitcativos.”. Ph.D. thesis, Universidade do Estado do Rio de Janeiro, Nova Friburgo, RJ (2018).
  • 16
    G.E. Sjoden . Deterministic adjoint transport applications for He-3 neutron detector design. Annals ofNuclear Energy , 29 (2002), 1055-1071.
  • 17
    R. Vargas & M.T. Vilhena. A closed-form solution for one-dimensional radiative condutive problem by the decomposition and LTSn methods. Journal of Quantitative Spectroscopy and Radiative Transfer , 61 (1999), 303-308.

Datas de Publicação

  • Publicação nesta coleção
    30 Nov 2020
  • Data do Fascículo
    Sep-Dec 2020

Histórico

  • Recebido
    26 Jul 2018
  • Aceito
    07 Jul 2020
Sociedade Brasileira de Matemática Aplicada e Computacional Rua Maestro João Seppe, nº. 900, 16º. andar - Sala 163 , 13561-120 São Carlos - SP, Tel. / Fax: (55 16) 3412-9752 - São Carlos - SP - Brazil
E-mail: sbmac@sbmac.org.br