Solução da Equação de Transporte de Nêutrons via Redes Neurais Informadas por Física: uma abordagem pedagógica

Solution of Neutron Transport Equation by Physics Informed Neural Network: a pedagogial approach

Resumo

Neste trabalho, trazemos um material pedagógico sobre a equação de transporte de nêutrons aplicada na análise de reatores nucleares. A apresentação é dividida em dois momentos. No primeiro, é ralizada uma revisão sobre alguns conceitos de física nuclear, com a finalidade de deduzir a equação de transporte de nêutrons, e assim, na sequência resolver analiticamente tal equação para o reator nuclear crítico com formato cilíndrico. No segundo momento, uma revisão sobre redes neurais informadas por física é apresentada. Posteriormente a equação de transporte de nêutrons para o reator nuclear com formato cilíndrico em funcionamento crítico é resolvida com o uso das PINNs. Por fim, a solução analítica é comparada com a solução obtida por intermédio das PINNs, evidenciando que o método de solução via redes neurais é de fato promissor.

Palavras-chave:
Inteligência artificial; aprendizado de máquina; redes neurais informadas por física; PINNs; equação de transporte de nêntrons

Abstract

In this work, we present pedagogical text on the neutron transport equation as applied to the analysis of nuclear reactors. The presentation is divided into two main parts. In the first part, we review fundamental concepts of nuclear physics with the aim of deriving the neutron transport equation. We then proceed to solve this equation analytically for critical nuclear reactors with cylindrical geometry. In the second part, we provide a review of Physics-Informed Neural Networks (PINNs). Subsequently, the neutron transport equation for the reactor configuration previously discussed is solved using the PINN approach. Finally, the analytical solution is compared with the solution obtained through the use of PINNs, demonstrating that the neural network-based method is indeed a promising approach.

Keywords:
Artificial inteligence; machine learning; Physics-Informed Neural Networks; PINNs; neutrons transport equation

1. Introdução

As tecnologias nucleares, apesar de serem bastante mitificadas, fazem parte do nosso cotidiano. Seja por meio de um exame de radiodiagnóstico, ou pelo consumo de alimentos conservados com uso de radiação, ou mesmo pela utilização de energia elétrica oriunda de uma usina nuclear, estamos imersos no universo da física nuclear. Dessa forma, compreedê-la nos ajuda a mitigar preconceitos existentes acerca da área nuclear e nos proteger contra o mau uso dessa fonte energética. Contudo, materiais instrucionais sobre energia nuclear escritos em língua portuguesa são escassos, sobretudo os destinados ao público leigo. Em particular, apesar de constar na literatura de língua inglesa bons materiais didáticos que discutem física de raetores nucleares, como por exemplos as referências [¹, ², ³, ⁴, ⁵, ⁶, ⁷, ⁸], tais bibliografias não encontram paralelo em língua portuguesa, comprometendo em parte, o aprendizado de alunos de graduação, bem como de interessados no assunto. Nessa perspectiva, este texto traz, num primeiro momento, uma revisão pedagógica acerca da equação de transporte de nêutrons, enfatizando os reatores nucleares críticos. Assim, soluções da equação de transporte nêutrons são apresentadas para o reator em formato cilíndrico. Tais soluções abarcam dois momentos: no primeiro momento são desenvolvidas analiticamente, e posteriormente são solucionadas mediante o uso de redes neurais informadas por física (PINNs).

O uso das PINNs (Physics-Informed Neural Networks) teve início de forma mais estruturada e reconhecida por volta de 2017, com destaque para os trabalhos pioneiros de Maziar Raissi, Paris Perdikaris e George Karniadakis. Esses pesquisadores propuseram formalmente a ideia de incorporar equações diferenciais (como leis físicas) diretamente no treinamento de redes neurais, de modo que a solução obtida respeitasse essas leis – mesmo em regiões sem dados. O trabalho de maior destaque nesse contexto, intitulado Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, foi escrito em 2019 pelo trio consituído por Raissi, Perdikaris e Karniadakis [⁹]. Apesar de existirem estudos anteriores que combinavam aprendizado de máquina e modelagem física, foi esse trabalho que consolidou o termo PINNs e abriu caminho para uma ampla adoção da técnica em áreas como fluidodinâmica, física quântica, engenharia biomédica, dentre outras [¹⁰, ¹¹, ¹²]. Desde então, o campo tem crescido rapidamente, com muitos avanços teóricos e práticos.

Nesse caminho, as Redes Neurais Informadas por Física representam um avanço significativo na interface entre aprendizado de máquina e modelagem científica. Sua importância reside na capacidade de incorporar diretamente leis físicas, expressas por equações diferenciais, dentro da estrutura de redes neurais, garantindo que as soluções obtidas respeitem os princípios fundamentais da física. Essa abordagem tem sido aplicada com sucesso em diversas áreas, sobretudo em física de reatores nucleares [¹³, ¹⁴, ¹⁵]. As PINNs são especialmente úteis para resolver equações diferenciais parciais e ordinárias em domínios complexos ou com dados escassos, onde métodos numéricos tradicionais enfrentam dificuldades. Ao conciliar dados empíricos e conhecimento físico, elas oferecem uma ferramenta robusta e flexível para simulações mais precisas e eficientes em problemas científicos e de engenharia. Nesse caminho, a nossa meta aqui foi incorporar, de uma forma simples e pedagógica, as PINNs na solução da equação de transporte de nêutrons para reatores nucleares críticos no formato cilíndrico; e no fim, comparar as soluções obtidas via PINNs com as soluções analíticas corriqueiras, também apresentadas neste trabalho. Com isso, almejamos que este trabalho auxilie o leitor, sobretudo o indivíduo lusófono, na compreensão de temas mais complexos tanto do arcabouço relativo à física de reatores nucleares, como no estudo de redes neurais informadas por física. Essa expectativa justifica-se em decorrência do trabalho apresentar uma sólida e cuidadosa revisão sobre a física de nêutrons, apresentando a equação de transporte aplicada no contexto de reatores nucleares, bem como por ter trazido uma discussão muito detalhada sobre PINNs, destinada principalmente a que está sendo introduzido no assunto.

Nessa perspectiva, a apresentação deste trabalho segue a sequência: na seção 2, uma revisão detalhada sobre a equação de transporte de nêutrons é destacada, trazendo ao leitor conceitos importantes como seção de choque e fluxo de nêutrons; na seção 3, aplicamos a equação de transporte de nêutrons na análise de reatores nucleares críticos de formato cilíndrico, obtendo, de forma detalhada, as soluções analíticas; na seção 4, apresentamos uma revisão sobre redes neurais informadas por física; na seção 5, trazemos as soluções para a equação de transporte obtidas com o uso de PINNs; por fim, na seção 6, as considerações finais e perspectivas são elencadas.

2. A Equação de Transporte de Nêutrons

Nesta seção, apresentaremos uma dedução da equação de transporte de nêutrons elencando algumas consequências físicas oiundas de sua interpretação.

A difusão de qualquer substância (neste texto, entenda-se por substância qualquer elemento que flua, seja um gás ou mesmo um grupo de partículas) é um problema já bem delineado pela física. Nesse sentido, o formalismo geral que será apresentado aqui seguirá de perto problemas análogos como a propagação do calor ou de um gás através de um dado meio. Nesse ideário, o leitor já deve ter vivenciado a situação de estar numa sala onde é aberto um recipiente de perfume. Nos instantes iniciais, após a abertura do frasco, apenas as pessoas que estão mais próximas sentirão o cheiro do perfume, mas após algum tempo, pessoas que estejam mais afastadas do recipiente poderão sentir o odor da substância. Isso ocorre porque o odor se difunde pela sala, e sua difusão obedece algumas leis e equações, como a lei de Fick e a equação de transporte, que veremos na sequência. O problema da difusão de nêutron é similar ao do perfume, ou seja, os nêutrons também se difundem por um dado meio, e o comportamento dessas partículas enquanto se difundem pode ser modelado por equações análogas à do perfume. Contudo, no caso dos nêutrons aparecem alguns ingredientes adicionais, pois a colisão de nêutrons com algumas partículas, como núcleos de urânio, podem levar ao surgimento de outros nêutrons, e essa particularidade deve ser incluída nas equações, o que será feito na dedução das equações próprias desse sistema.

2.1. O fluxo de nêutrons

A equação que descreve a difusão de uma dada substância representa a evolução temporal que descreve alguma característica da substância. Por exemplo, a equação de difusão do calor é representada pela evolução da quantidade de calor num dado ponto do espaço num dado instante. No caso da difusão dos nêutrons, a grandeza física que será utilizada como função a ser determinada é o fluxo de nêutrons, que será descrito a seguir.

O fluxo de nêutrons $\varphi$ é definido por

(1) $$\varphi =n\cdot v,$$

em que $n$ representa a densidade de nêutrons, isto é, o número de nêutrons por centímetro cúbico numa dada região, e $v$ representa a velocidade média dos nêutrons (consideraremos neste rabalho apenas nêutrons monoenergéticos, isto é, com a mesma velocidade). Perceba que a unidade do fluxo de nêutrons é nêutrons por centímetro quadrado por segundo ($\text{nêutrons}/c{m}^{2}s$). O fluxo de nêutrons, em geral, é uma grandeza que depende da posição espacial e do tempo.

Em física nuclear, é usual trabalhar com uma grandeza denominada seção de choque, denotada por $\sigma$, a qual fornece a probabilidade de uma partícula colidir com a outra. A seção de choque é geralmente determinada experimentalmente, e depende de alguns parâmetros, como por exemplo da energia cinética (ou velocidade) da partícula incidente (é usual adotar um referencial em que uma das partículas esteja em movimento (a partícula incidente) e a outra esteja parada (partícula alvo). Assim, em materiais de física nuclear, há tabelas que contêm seções de choques para vários pares de partículas com diferentes velocidades. Na essência, a seção de choque (ou seção eficaz, em inglês cross section) representa a probabilidade de ocorrência de uma interação entre um nêutron e um núcleo atômico. Essa interação pode ser:

• Dispersão elástica ou inelástica (o nêutron colide e desvia, com ou sem perda de energia),

• Absorção (o nêutron é capturado),

• Fissão (em núcleos fissionáveis como o U-235).

A seção de choque é medida em unidades de área, sendo o barn ($b$) a unidade padrão:

$$1\text{barn}={10}^{-24}c{m}^2.$$

A seção de choque total é a soma de todas as seções de choque parciais:

$${\sigma }_t={\sigma }_s+{\sigma }_a,$$

em que ${\sigma }_t$ representa a seção de choque total, ${\sigma }_s$ representa a seção de choque de dispersão ou espalhamento (scattering), e ${\sigma }_a$ representa a seção de choque de absorção (inclui fissão, captura, etc.). Vale a pena enfatizar que a seção de choque de fissão é uma subcategoria da seção de choque de absorção, ou seja, a seção de choque de fissão é a soma das seções de choque de fissão, captura radiativa, dentre outras.

A seguir, na Tabela 1, apresentamos alguns valores das seções de choque dos isótopos U-235 (urânio-235) e U-238 (urânio-238) para nêutrons térmicos (ou lentos), isto é, nêutrons com energia cinética de $\mathrm{0,025}eV$.

Tabela 1
Valores aproximados das seções de choque dos isótopos U-235 e U-238. Fonte: [⁸].

Isótopo	Tipo de Interação	$\sigma$ (barns)
U-235	Fissão	$\approx \mathrm{582,6}$
	Captura Radiativa	$\approx \mathrm{98,8}$
	Absorção Total	$\approx \mathrm{681,4}$
	Espalhamento	$\approx 14$
	Total	$\approx \mathrm{695,4}$
U-238	Fissão	$\approx \mathrm{0,002}$
	Captura Radiativa	$\approx \mathrm{2,68}$
	Absorção Total	$\approx \mathrm{2,682}$
	Espalhamento	$\approx \mathrm{8,871}$
	Total	$\approx \mathrm{11,55}$

Na Tabela 1, como pode ser percebido, a seção de choque de absorção total é a soma seção de choque de fissão com a seção de choque de captura radiativa, enquanto a seção de choque total é igual a soma da seção de choque de absorção total com a seção de choque de espalhamento. O U-235 é fissionável com nêutrons térmicos, por isso tem uma seção de choque de fissão muito alta a baixas energias. Já o U-238 só sofre fissão com nêutrons rápidos, e sua seção de choque de fissão térmica é praticamente desprezível. Esses valores mudam bastante em energias mais altas (região epitérmica ou rápida).

A Figura 1 traz o gráfico das seções de choque dos isótopos U-235 e U-238, mostrando a dependência desse parâmetro com a temperatura ou com a energia cinética dos nêutrons.

Figura 1
Dependência da seção de choque dos isótopos U-235 e U-238 em função da energia cinética dos nêutrons (Fonte: autoria própria).

No gráfico, as linhas contínuas mostram a seção de choque total, enquanto as linhas tracejadas mostram a seção de choque de fissão. Conforme podemos notar, o U-235 tem uma seção de fissão muito alta em nêutrons térmicos (à esquerda do gráfico), o que o torna ideal para reatores nucleares térmicos. Enquanto isso, o U-238 praticamente não sofre fissão com nêutrons térmicos, mas a seção de fissão aumenta com a energia, tornando-o relevante em reatores rápidos.

Se multiplicarmos a densidade atômica $N$ pela seção de choque $\sigma$, obtemos um parâmetros denominado seção de choque macroscópica, denotada por $\text{Σ}$,

(2) $$\text{Σ}=N\sigma .$$

Assim como a seção de choque $\sigma$ é definida para cada tipo de interação, a seção de choque macroscópica $\text{Σ}$ também. Por exemplo, a seção de choque de absorção do U-235 e U-238 para nêutrons com energia cinética de $\mathrm{0,025}eV$ são $\mathrm{681,4}b$ e $\mathrm{2,682}b$, respectivamente. Como as densidades atômicas de U-235 e U-238 no urânio natural são $\mathrm{3,48}\cdot {10}^{\text{20}}\text{átomos}/c{m}^3$ e $\mathrm{0,0483}\cdot {10}^{\text{24}}\text{átomos}/c{m}^3$, respectivamente; então, a seção de choque macroscópica de absorção do U-235 e U-238 para os nêutrons térmicos são, respectivamente, $\mathrm{0,237}c{m}^{-1}$ e $\mathrm{0,130}c{m}^{-1}$. Assim, a seção de choque macroscópica de absorção do urânio natural é dada por $\mathrm{0,367}c{m}^{-1}$.

Perceba que se conhecemos o fluxo de nêutrons, isto é, o número de nêutrons em cada centímeto quadrado de uma dada região em cada segundo, e também a seção de choque total das partículas-alvo disponíveis nesta região, podemos calcular a quantidade de colisões que ocorrerão por centímetro cúbico por segundo nesta região mediante o produto

$$F={\text{Σ}}_t\varphi ,$$

em que $F$ é denominado densidade de colisões. Esse é mais um parâmetro relevante na engenharia nuclear. De posse desses parâmetros, podemos discutir o próximo ponto: a lei de Fick.

2.2. A Lei de Fick

A teoria da difusão é baseada na lei de Fick. Tal lei foi originalmente utilizada na difusão química. É de conhecimento que quando, numa solução, a concentração de determinado soluto numa região é maior do que em outra, o soluto se difunde da região de maior concentração para a de menor concentração. Além disso, sabe-se também que a taxa com que o soluto flui é proporcional ao negativo do gradiente da concentração do soluto. Os nêutrons, em certas condições, possuem o comportamento similar ao do soluto numa solução. Isto é, se a desnidade ou o fluxo de nêutrons for maior numa região de um reator do que em outra, existirá uma corrente de nêutrons fluindo da região de maior densidade para a de menor densidade. Assim, se considerarmos que o fluxo varia ao longo da direção $x$, a lei de Fick fica dada por [⁶, ⁸]

(3) $$J_x=-D\frac{d\varphi }{dx},$$

em que $D$ é o coeficiente de difusão, o qual depende de parâmetros relacionados ao sistema analisado; no caso considerado, $D$ depende dos parâmetros associados ao meio em que os nêutrons estão fluindo. $J_x$ é denominada densidade de corrente de nêutrons, sendo uma medida do número de nêutrons que passa, por unidade de tempo, através de uma unidade de área perpendicular à direção $x$. A unidade de $J_x$ é a mesma do fluxo, isto é, $\text{nêutrons}/c{m}^{2}s$.

Se o fluxo é função das três coordenadas espaciais, definimos

(4) $$\overrightarrow{J}=-D\nabla \varphi ,$$

em que $\overrightarrow{J}$ é denominado vetor densidade de corrente de nêutrons. O significado físico de $\overrightarrow{J}$ pode ser obtido se considerarmos o produto escalar

$$\overrightarrow{J}\cdot \widehat{n}=J_n,$$

onde $\widehat{n}$ é o vetor normal unitário numa direção arbitrária. $J_n$ obtido é igual à quantidade de nêutrons por segundo que atravessam uma unidade de área normal à direção definida por $\widehat{n}$.

A fim de compreendermos melhor esses conceitos, consideremos que uma fonte de nêutrons pontual que emita $Q$ nêuntrons por segundo possua o fluxo de nêutrons dado por

$$\varphi \left(r\right)=\frac{Q{e}^{-r/L}}{4\pi Dr},$$

em que $r$ é a distância radial da fonte até algum ponto, $D$ é o coeficiente de difusão, e $L$ é uma constante. A partir deste fluxo, podemos calcular a corrente. Notemos que $\varphi$ depende apenas de $r$. Assim, se tomarmos o operador gradiente em coordenadas esféricas, considerando apenas a componente radial, temos que $\nabla =\widehat{r}\frac{d}{dr}$, em que $r$ é o vetor unitário na direção radial. Assim, a densidade de corrente, que também só terá componente na direção radial, é dada por

$$\overrightarrow{J}=-D\nabla \varphi ,$$

que no caso em tela é dado por

$$\overrightarrow{J}=-D\widehat{r}\frac{d}{dr}\left(\frac{Q{e}^{-r/L}}{4\pi Dr}\right),$$

a partir de que obtemos

$$\overrightarrow{J}=\widehat{r}\frac{Q}{4\pi }\left(\frac{1}{r^2}+\frac{1}{rL}\right)e^{-r/L}.$$

A partir de $J$, podemos calcular a quantidade de nêutrons por segundo que atravessam uma unidade de área normal à direção definida por $\widehat{r}$. Para isso, notemos que no caso considerado $\overrightarrow{J}$ é sempre perpendicular à superfície de uma esfera. Assim, $\overrightarrow{J}\cdot \widehat{r}$ nos fornece a quantidade de nêutrons que atravessam cada unidade de área da esfera. E assim, a quantidade de nêutrons que, a cada segundo, atravessam a esfera inteira é dada por $4\pi r^2\overrightarrow{J}\cdot \widehat{r}$, o que nos fornece $Q\left(1+\frac{r}{L}\right)e^{-r/L}$.

Para finalizar esta subseção, vamos discutir o significado do coeficiente $D$. Lembrando que $\overrightarrow{J}$ e $\varphi$ possuem a mesma unidade e visualizando a Eq. (3), percebemos que $D$ tem que ter dimensão de comprimento. Pode-se mostrar (o que é muito longo para ser feito aqui neste texto introdutório), que o coeficiente de difusão pode ser escrito como [⁸]

(5) $$D=\frac{1}{3{\text{Σ}}_S\left(1-\overline{\mu }\right)},$$

em que $\overline{\mu }$ é o valor médio do cosseno do ângulo que os nêutrons são espalhados após a colisão com as partículas-alvo. Para os núcleos mais usados nos reatores, $\overline{\mu }$ pode ser calculado por meio de

$$\overline{\mu }=\frac{2}{3A},$$

onde $A$ é o número de massa atômica da partícula-alvo. Se definirmos o livre caminho médio de transporte ${\lambda }_{tr}$ como

$${\lambda }_{tr}=\frac{1}{{\text{Σ}}_S\left(1-\overline{\mu }\right)},$$

podemos reescrever $D$ na forma simples

$$D=\frac{1}{3{\lambda }_{tr}}.$$

2.3. A equação de continuidade

Outro ingrediente fundamental para a obtenção da equação de transporte de nêutrons é a equação de continuidade. Essa equação aparece nas mais variadas áreas da física, de forma que o leitor já possa ter se deparado com ela em outros contextos. A fim de obtê-la para o caso dos nêutrons, consideremos um volume arbitrário $V$ no qual há nêutrons e outras partículas. Com o passar do tempo, a quantidade de nêutrons dentro de $V$ pode mudar por algum dos seguintes processos: se nêutrons escaparem para fora de $V$, se nêutrons forem absorvidos por alguma pártícula presente em $V$, ou se fontes emissoras de nêutrons estão presentes em $V$. A equação de continuidade, a qual será obtida nesta subseção, é nada mais que a linguagem matemática para dizer que, ao levarmos em conta os três processos descritos, a quantidade de nêutrons deve se conservar.

Resumidamente, podemos escrever que a taxa de variação do número de nêutrons em V é igual à taxa com que os nêutrons são produzidos em $V$, subtraído da taxa que os nêutrons são absorvidos em $V$, subtraído também da taxa com que os nêutrons escapam de $V$. Temos que contabilizar cada termo sugerido neste resumo.

Com a finalidade de contabilizarmos a variação do número de nêutrons com o tempo, consideremos $n$ como a densidade de nêutrons em qualquer ponto em $V$. Assim, o total de nêutrons em $V$ é dado por

$${\int }_{V}n(\overrightarrow{r},t)\mathit{d}V,$$

onde a integração é realizada em todo o volume $V$.

A taxa de variação do número de nêutrons é dada por

$$\frac{d}{dt}{\displaystyle \underset{V}{\int }ndV}={\displaystyle \int \frac{\partial n}{\partial t}dV.}$$

Para computarmos a taxa de produção, consideremos que os nêutrons são emitidos por uma fonte localizada em $V$ que produza $s$ nêutrons por centímeto cúbico a cada segundo. Assim, a taxa de produção é dada por

$${\displaystyle \underset{V}{\int }sdV.}$$

Com essa integral obtemos o número de nêutrons produzidos no volume $V$ a cada segundo.

A taxa de absorção por centímetro cúbico por segundo é dada por ${\text{Σ}}_a\varphi$, em que ${\text{Σ}}_a$ é a seção de choque macroscópica da partícula-alvo (essa partícula-alvo, que no reator nuclear pode ser o núcleo de U-235, poderá absorver o nêutron após a colisão). Assim, a taxa de absorção é dada por

$$\underset{V}{{\displaystyle \int }}{\Sigma }_a\varphi dV.$$

Após o cálculo dessa integral, obtemos o número de nêutrons absorvidos dentro de $V$ a cada segundo.

A taxa de escape de nêutrons é encontrada por meio da densidade de corrente de nêutrons. Para isso, considere $\overrightarrow{J}$ como o vetor densidade de corrente na superfície de $V$ e $\widehat{n}$ é o vetor normal unitário direcionado para fora da superfície. Assim, sabemos que $\overrightarrow{J}\cdot \widehat{n}$ nos fornece o número de nêutrons que atravessam a superfície $A$ que delimita $V$ - de dentro para fora - por $c{m}^2/s$. Com isso, a taxa de vazamento é dada por

$$\underset{A}{{\displaystyle \int }}\overrightarrow{J}\cdot \widehat{n}dA.$$

Como $A$ é uma superfície fechada, podemos usar o teorema da divergência e reescrever a taxa de vazamento como

$$\underset{V}{{\displaystyle \int }}\nabla \cdot \overrightarrow{J}dV.$$

Nesse caminho, a taxa com que a quantidade de nêutrons dentro de $V$ varia com o tempo é dada por

$$\underset{V}{{\displaystyle \int }}\frac{\partial n}{\partial t}dV=\underset{V}{{\displaystyle \int }}sdV-\underset{V}{{\displaystyle \int }}{\Sigma }_a\varphi dV-\underset{V}{{\displaystyle \int }}\nabla \cdot \overrightarrow{J}dV.$$

Como todas as integrais estão relacionadas ao mesmo volume, podemos escrever

(6) $$\frac{\partial n}{\partial t}=s-{\text{Σ}}_a\varphi -\nabla \cdot \overrightarrow{J},$$

que é a forma geral da equação de continuidade. Nos casos em que a densidade de nêutrons não é função do tempo, obtemos

(7) $$\nabla \cdot \overrightarrow{J}+{\text{Σ}}_a\varphi -s=0,$$

que é equação de continuidade do estado estacionário.

2.4. A equação de transporte

Se substituirmos a Eq. (4) na Eq. (6), obtemos

(8) $$D{\nabla }^2\varphi -{\text{Σ}}_a\varphi +s=\frac{\partial n}{\partial t},$$

que é a equação de transporte de nêutrons.

Para o caso estacionário, obtemos

(9) $$D{\nabla }^2\varphi -{\text{Σ}}_a\varphi =-s,$$

a qual pode ser escrita na forma

(10) $${\nabla }^2\varphi -\frac{1}{L^2}\varphi =-\frac{s}{D},$$

em que $L=\sqrt{D/{\text{Σ}}_a}$ é denominado comprimento de difusão, termo que aparece com frequência na engenharia nuclear.

A Eq. (10) é a equação de difusão de nêutrons estacionária mais genérica, ainda não adaptada para a análise de um retator nuclear. Em breve faremos essa adaptação. Antes disso, vamos aplicá-la a uma situação prática.

2.4.1. Uma placa infinita irradiando nêutrons

Consideremos uma fonte de nêutrons que seja um plano infinito que emite $Q_0$ nêutrons por $c{m}^2/s$. Se considerarmos que o plano é o $y-z$, localizado em $x=0$, podemos admitir que o fluxo é uma função apenas da coordenada $x$, que é a distância de onde o fluxo é medido até o plano. Assim, podemos considerar que não há fontes em nenhum outro ponto além de $x=0$, ou seja, utilizaremos a fonte no plano como uma condição de contorno. Dessa forma, a Eq. (10) pode ser escrita para este caso como

(11) $${\nabla }^2\varphi -\frac{1}{L^2}\varphi =0,x\ne 0.$$

A solução geral da Eq. (11) é dada por

$$\varphi \left(x\right)=A{e}^{-x/L}+B{e}^{x/L},$$

em que $A$ e $B$ são constantes a serem determinadas com o uso das condições de contorno. Por ora, consideraremos apenas $x>0$. Assim, para que o fluxo seja finito, o coeficiente $B$ deve ser nulo. Dessa forma, o fluxo de nêutrons fica dado por

$$\varphi \left(x\right)=A{e}^{-x/L}.$$

O coeficiente $A$ é determinado com o uso do fato de que na origem, em $x=0$, há uma fonte que emite $Q_0$ nêutrons por centímetro quadrado por segundo. Como o plano emite nos dois sentidos, admitamos, por simetria, que no sentido $x>0$ a emisão seja de $Q_0/2$ nêutrons por centímetro quadrado por segundo. Dessa forma, a condição de contorno nos diz que

$$\underset{x\to 0}{lim}J\left(x\right)=\frac{Q_0}{2}.$$

A partir da lei de Fick, podemos escrever

$$J\left(x\right)=-D\frac{d\varphi }{dx}=\frac{DA}{L}{e}^{-x/L}.$$

Logo, usando a condição de contorno, temos que

$$A=\frac{Q_{0}L}{2D}.$$

E assim, o fluxo para $x>0$ é dado por

$$\varphi \left(x\right)=\frac{Q_{0}L}{2D}{e}^{-x/L}.$$

Se quisermos uma solução que valha também para $x<0$, por simetria, basta que escrevamos

(12) $$\varphi \left(x\right)=\frac{Q_{0}L}{2D}{e}^{-\left|x\right|/L}.$$

Os resultados deste exemplo serão utilizados num caso mais complexo que trataremos a seguir.

2.4.2. Uma placa infinita irradiando nêutronsnuma região fechada

No próximo exemplo, consideraremos o plano infinito do exemplo anterior, irradiando nêutrons a uma taxa $Q_0$ nêutrons por centímetro quadrado por segundo, envolto pelos dois lados por paredes infinitas. Consideremos que o plano que irradia esteja localizado na origem do sistema de coordenadas, em $x=0$, e a uma ditância $d$ de ambas as paredes. Ou seja, a parede da esquerda está localizada em $x=-d$, enquanto a da direita está em $x=d$. A equação que devemos resolver é a mesma do exemplo anterior, Eq. (11) pois consideraremos a placa irradiando como uma condição de contorno, isto é,

$$\underset{x\to 0}{lim}J\left(x\right)=\frac{Q_0}{2}.$$

A solução é dada por $\varphi \left(x\right)=A{e}^{-x/L}+B{e}^{x/L}$, e aqui não podemos assumir nenhuma das duas constantes como nula, já que o sistema não inclui o infinito, pois é limitado pelas paredes.

Agora, além dessa condição, temos também que o fluxo de nêutrons nas paredes deve ser nulo, ou seja,

$$\varphi \left(d\right)=\varphi \left(-d\right)=0.$$

Assim, temos então

$$\varphi \left(d\right)=A{e}^{-d/L}+B{e}^{d/L}=0,$$

e

$$\varphi \left(-d\right)=A{e}^{d/L}+B{e}^{-d/L}=0,$$

Com isso, podemos escrever $B=-A{e}^{-2d/L}$, e a solução será

$$\varphi \left(x\right)=A\left[e^{-x/L}-e^{\left(x-2d\right)/L}\right].$$

A constante $A$ pode ser determinada, como no exemplo anterior, com a condição da fonte $Q_0$. Usando a lei de Fick obtemos:

$$J\left(x\right)=\frac{DA}{L}\left[e^{-x/L}+e^{\left(x-2d\right)/L}\right].$$

E assim,

$$A=\frac{Q_{0}L}{2D}{\left[1+e^{-2d/L}\right]}^{-1}.$$

Logo, para $x>0$, obtemos

$$\varphi \left(x\right)=\frac{Q_{0}L}{2D}\left[\frac{e^{-x/L}-e^{\left(x-2d\right)/L}}{1+e^{-2d/L}}\right].$$

Por simetria, a solução válida para qualquer valor de $x$ é dada por

(13) $$\varphi \left(x\right)=\frac{Q_{0}L}{2D}\left[\frac{e^{-\left|x\right|/L}-e^{\left(\left|x\right|-2d\right)/L}}{1+e^{-2d/L}}\right].$$

Multiplicando o numerador e o denominador da Eq. (13) por $e^{d/L}$, chegamos a

(14) $$\varphi \left(x\right)=\frac{Q_{0}L}{2D}\left[\frac{e^{-\left(\left|x\right|-d\right)/L}-e^{\left(\left|x\right|-d\right)/L}}{e^{d/L}+e^{-d/L}}\right].$$

Agora, usando as funções trigonométricas hiperbólicas, isto é $2\cosh w=e^w+e^{-w}$ e $2\sinh w=e^w-e^{-w}$

, a Eq. (14) passa a ser escrita na forma

(15) $$\varphi \left(x\right)=\frac{Q_{0}L}{2D}\left[\frac{\sinh \left[\left(d-\left|x\right|\right)/L\right]}{\cosh \left(d/L\right)}\right].$$

Neste ponto, é importante comentarmos que a condição de contorno $\varphi \left(d\right)=\varphi \left(-d\right)=0$ não é uma condição física, pois o fluxo de nêutrons não é realmente nulo nas paredes; caso fossem, não haveria corrente de nêutrons através das paredes, e nenhum nêutron seria encontrado além delas. Assim, na verdade $\varphi \approx 0$ nas paredes. Sob um ponto de vista mais formal, a diferença reside na condição de fronteira que adotamos, ou seja, utilizamos a condição de Dirichlet, na qual o nó se encontra rigorosamente na fornteira. A abordagem fica fisicamente mais aceitável caso utilizássemos a condição de Marshak [¹⁶], na qual o nó encontra-se fora do domínio físico, isto é, numa fronteira extrapolada. De forma simplificada, na condição de Dirichlet, o fluxo cai abruptamente para zero na parede, enquanto na de Marshak, o fluxo atinge zero um pouco depois da parede. Em síntese, na condição de Marshak, a extrapolação da fornteira pode ser escrita como $\varphi \left(d\pm \delta \right)=0$, em que $\delta \approx 2D$. O valor mais preciso da teoria de transporte é $\delta =\mathrm{0,7104}\cdot \left(3D\right)\approx \mathrm{2,13}D$. Esse resultado advém da condição

$$\varphi +2D\frac{d\varphi }{dx}=0,$$

o que implica também $J=\varphi /2$. Optamos por não deduzir aqui os detalhes do último resultado mostrado. Na prática, continuaremos usando a condição de Dirichlet, mas compreendendo que, fisicamente, há fluxo de nêutron na parede, o que possibilita a corrente de nêutrons através dela, e assim, o vazamento de nêutrons nos casos reais. Assim, com o intuito de simplificar os cálculos, continuaremos a utilizar a condição de Dirichlet.

Com tal explicação concluída, podemos retomar a discussão prática sobre os planos infinitos. A Eq. (15) nos fornece bastante informação acerca do sistema considerado. Por exemplo, suponha que queiramos determinar o número de nêutrons que escapa do sistema por segundo de cada $c{m}^2$ através das duas paredes. Para isso, é suficiente calcularmos a corrente de nêutrons na parede, ou seja, em $x=d$ e em $x=-d$ e adicionarmos os dois resultados.

Fazendo isso, a partir da Eq. (15) podemos deduzir que

$$J\left(x\right)=-\frac{Q_0\cosh \left[\left(d-\left|x\right|\right)/L\right]}{2\cosh \left(d/L\right)}.$$

Logo,

$$J\left(d\right)=-\frac{Q_0}{2\cosh \left(d/L\right)}.$$

Da mesma forma,

$$J\left(-d\right)=\frac{Q_0}{2\cosh \left(d/L\right)}.$$

E assim, a quantidade de nêutrons que escapa do sistema por segundo de cada $c{m}^2$ através das duas paredes é igual a $J\left(d\right)+J\left(-d\right)$, ou seja,

$$J\left(d\right)-J\left(-d\right)=\frac{Q}{\cosh \left(d/L\right)}.$$

A partir deste resultado, podemos determinar a probabilidade com que os nêutrons emitidos possam vazar pelas paredes. Se a fonte emite $Q$ nêutrons por centímetro quadrado por segundo e se escapam $\frac{Q_0}{\cosh \left(d/L\right)}$, essa probabilidade é dada por

$$\frac{J\left(d\right)+J\left(-d\right)}{Q}=\frac{1}{\cosh \left(d/L\right)}.$$

Esses cálculos, apesar de serem relativamente simples, são muito importantes na elaboração de projetos de blindagem, quando se quer reduzir a quantidade de nêutrons que escapam de determinado sistema.

2.5. A equação de transporte de nêutronsno contexto do reator nuclear

Nesta subseção, deduziremos a equação de transporte de nêutrons apropriada para se analizar o fluxo de nêutrons num reator nuclear. Consideraremos o modelo mais simples que simula um reator nuclear: um reator rápido e crítico contendo uma mistura de combustível e refrigerante, consistindo de apenas uma região. Este reator é denominado reator nú e a dinâmica dos nêutrons em seu interior podem ser descritos pela Eq. (10), isto é,

$$D{\nabla }^2\varphi -{\text{Σ}}_a\varphi =-s,$$

No reator crítico, o que estamos considerando aqui, são os nêutrons emitidos na fissão. Aqui vale recordar que quando um nêutron atinge um núcleo de urânio-235, ele pode ser absorvido, formando um núcleo instável de urânio-236, que se divide em dois fragmentos menores (núcleos fi-lhos) e libera mais nêutrons, além de uma grande quantidade de energia.

A reação mais geral pode ser escrita como

$${}^{235}U{+}^{1}n\to {}_{}{}^{236}U_{}^{\ast }\to \text{Fragmentos da Fissão}+{\left(\text{2 ou 3}\right)}^{1}n+\text{Energia}.$$

Um exemplo específico de fissão é

$$\begin{array}{l}{}^{235}U{+}^{1}n{\to }^{236}{U}^{\ast }\to \\ {}^{92}Kr{+}^{141}Ba+3^{\text{1}}n+Energia.\end{array}$$

A energia liberada é em torno de $200MeV$ por fissão. Dessa forma, a fonte de nêutrons $s$ pode ser dada por

$$s=\eta {\text{Σ}}_{aF}\varphi ,$$

em que $\eta$ é o número médio de nêutrons emitidos por fissão, e ${\text{Σ}}_{aF}$ é a seção macroscópica de choque total do combustível.

Notemos que o termo de fonte pode ser reescrito como

$$s=\eta f{\text{Σ}}_a\varphi ,$$

em que a razão $f=\frac{{\text{Σ}}_{aF}}{{\text{Σ}}_a}$ é denominado fator de utilização e ${\text{Σ}}_a$ é a seção macroscópica de choque da mistura combustível e refrigerante. Assim, $f$ é igual a fração de nêutrons absorvidos no reator e que de fato são absorvidos no combustível. Ou seja, $f$ fornece a fração de nêutrons que de fato podem provocar a fissão no urânio-235.

O produto $\eta f$ é denominado fator de multiplicação de um meio infinito, representado por $k_{\infty }=\eta f$. Um meio infinito é definido como um meio em que não há perda de nêutrons por fuga.

De forma geral, o fator de multiplicação, geralmente representado por

$$k=\frac{\text{número de nêutrons na geração atual}}{\text{número de nêutrons na geração anterior}}.$$

$k$ é um parâmetro fundamental na física de reatores nucleares. Ele indica a razão entre o número de nêutrons na geração atual e o número de nêutrons na geração anterior. Esse parâmetro pode ser classificado em dois tipos, o fator de multiplicação efetivo, $k_{eff}$, que leva em conta todas as perdas de nêutrons no reator, como nêutrons que escapam do reator e nêutrons absorvidos sem causar fissão; e o fator de multiplicação infinito (já mencionado anteriormente), $k_{\infty }$.

O valor de $k$ pode ser interpretado com os seguintes critérios:

• $k<1$ – Subcrítico: A reação de fissão está diminuindo; o reator "morre"com o tempo.

• $k=1$ -– Crítico: A reação está estável. Cada geração de nêutrons mantém a próxima. É o estado desejado para operação contínua de um reator.

• $k>1$ -– Supercrítico: A reação está se intensificando. O número de nêutrons (e a potência) cresce com o tempo.

Com isso, o termo de fonte pode finalmente ser escrito como

$$s=k_{\infty }{\text{Σ}}_a\varphi .$$

Subsituindo na Eq. (10), obtemos

(16) $${\nabla }^2\varphi +B^2\varphi =0,$$

em que

(17) $$B^2=\frac{k_{\infty }-1}{L^2}.$$

Nas (equações 16) e (17) temos o termo $B$ que é referido em inglês pela palavra Buckling medida em cm^-2. Essas equações serão aplicadas na próxima seção para o estudo do reator nuclear cilíndrico.

3. Solução Analítica da Equação de Transporte de Nêutrons aplicadas a reatores nucleares

Nesta seção, o nosso objetivo é aplicar a Eq. (16) para analisar reatores nucleares. Consideraremos nesta abordagem analítica o reator no formato cilíndrico, que é o caso mais comum em usinas nucleares, sobretudo nas usinas PWR [⁵]. Esperamos que esta seção forneça ao leitor subsídios que o auxilie na compreensão da física de reatores nucleares, sobretudo dos conceitos envolvidos no projeto de retatores.

Nesse caminho, vamos considerar o reator cilíndrico de raio $a$ e altura $H$. Devido à simetria do problema, optamos por utilizar coordenadas cilíndricas $(r,\theta ,z)$. Em nossa análise, consideraremos o caso em que o fluxo de nêutrons não dependa da coordenada $\theta$. Assim, o laplaciano apenas com a dependência em $r$ e $z$ é dado por

$${\nabla }^2=\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial }{\partial r}\right)+\frac{{\partial }^2}{\partial z^2}.$$

Com isso, a Eq. (16) fica escrita na forma

(18) $$\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \varphi }{\partial r}\right)+\frac{{\partial }^2\varphi }{\partial z^2}+B^2\varphi =0.$$

Assim, como foi feito na seção precedente, vamos supor que a solução seja dada por $\varphi (r,z)=R\left(r\right)Z\left(z\right)$. Assim, obtemos as seguintes equações diferenciais ordinárias

(19) $$\frac{1}{r}\frac{d}{dr}\left(r\frac{dR}{dr}\right)+B_r^{2}R=0,$$

(20) $$\frac{d^{2}Z}{d{z}^2}=-B_z^{2}Z,$$

onde $B_r^2$ e $B_z^2$ são as constantes de separação, e satisfazem

(21) $$B^2=B_r^2+B_z^2.$$

A solução da Eq. (20) é dada por

$$Z\left(z\right)=A_1\cos \left(B_{z}z\right)+C_1\sin \left(B_{z}z\right).$$

Utilizando a condição de contorno que o fluxo deve ser nulo nas tampas do reator (bases do cilindro), $z=H/2$ e $z=-H/2$, encontramos que a solução deve ser uma função par que satisfaça $Z\left(H/2\right)=0$.

Isso nos fornece para o reator crítico

$$\cos \left(B_z\frac{H}{2}\right)=0,$$

ou seja,

(22) $$B_z=\frac{\pi }{H}.$$

E assim, a solução da parte em $z$ fica dada por

(23) $$Z\left(z\right)=A_1\cos \left(\frac{\pi }{H}z\right).$$

Agora retomamos à solução da parte radial. Após manipulação, a Eq. (19) pode ser reescrita na forma

(24) $$r^2\frac{d^{2}R}{d{r}^2}+r\frac{dR}{dr}+B_r^2{r}^{2}R=0$$

Essa equação é um caso particular da equação diferencial de Bessel, que é dada dada por [¹⁷]

(25) $$r^2\frac{d^{2}R}{d{r}^2}+r\frac{dR}{dr}+\left(B_r^2{r}^2-p^2\right)R=0.$$

Ou seja, a Eq. (24) é a equação diferencial de Bessel, Eq. (25), com $p=0$. Sua solução é bem conhecida na literatura, pois é dada em termos de [¹⁷]

(26) $$R\left(r\right)=A_2{J}_0\left(B_{r}r\right)+C_2{Y}_0\left(B_{r}r\right),$$

em que $A$ e $C$ são constantes, enquanto $J_0\left(B_{r}r\right)$ e $Y_0\left(B_{r}r\right)$ são, respectivamente, a função de Bessel e a função de Neumann de ordem zero. Todas as identidades e propriedades sobre as funções de Bessel que utilizamos podem ser encontradas na referência. Enquanto a função de Bessel $J_0\left(B_{r}r\right)$ é bem comportada para qualquer valor da variável, a função de Neumann $N_0\left(B_{r}r\right)$ diverge na origem, ou seja, $\underset{r\to 0^{+}}{lim}{Y}_0\left(B_{r}r\right)=\infty$. Como o nosso problema envolve a origem, ou seja, o núcleo do reator cilíndrico, localizado em $r=0$, faz parte do sistema analisado, o coeficiente $C_2$ da Eq. (26) deve ser igual a zero. Por isso, a solução da parte radial fica dadapor

(27) $$R\left(r\right)=A_2{J}_0\left(B_{r}r\right).$$

Usando a condição de contorno que o fluxo deve ser nulo na lateral do cilindro, em $r=a$, temos

$$J_0\left(B_{r}a\right)=0.$$

Como estamos interessados no reator crítico, apenas a primeira raiz da função de Bessel de ordem zero é relevante. Essa raiz é igual a $\mathrm{2,405}$ [¹⁷]. Assim,

(28) $$B_r=\frac{\mathrm{2,405}}{a}$$

E a solução da parte radial toma a forma

(29) $$R\left(r\right)=A_2{J}_0\left(\frac{\mathrm{2,405}r}{a}\right).$$

Com isso, a solução do problema fica a princípio dada por

(30) $$\varphi (r,z)=A{J}_0\left(\frac{\mathrm{2,405}r}{a}\right)cos\left(\frac{\pi z}{H}\right),$$

em que $A=A_1{A}_2$ é uma constante que determinaremos a seguir.

Assim, como comentamos no caso dos planos infinitos, a condição de Dirichlet, $\varphi =0$ na fronteira da região considerada (tampas e laterais do cilindro), é apenas uma condição matemática. Na prática, o fluxo de nêutrons não é nulo na fronteira, pois se assim fosse, nêutrons não vazariam do reator. Na verdade, temos que $\varphi \approx 0$, o que advém da condição de Marshak [¹⁶], isto é,

$$\varphi (R,z)+2D\frac{\partial \varphi (R,z)}{\partial r}=0,$$

nas laterais do cilindro ($r=R$), e

$$\varphi (r,\pm H/2)+2D\frac{\partial \varphi (r,\pm H/2)}{\partial z}=0,$$

nas tampas do cilindro ($z=\pm H/2$). Na prática, trocamos $R\to R+\delta$ e $H=H+2\delta$, em que $\delta =2D$, nas equações oriundas da solução do reator, inclusive na Eq. (17). Como este trabalho tem cunho meramente pedagógico, continuaremos a admitir a condição de Dirichlet (entendendo-a apenas como formalidade matemática), mas sabendo que na realidade o fluxo de nêutrons não se anula nas fronteiras da região analisada.

A potência, grandeza física que se conhece (a priori) ao se projetar um reator, pode ser calculada por meio da integral [⁸]

$$P=E_R{\text{Σ}}_f{\int }_V\varphi (r,z)\mathit{d}V,$$

o que no caso considerado é igual a

(31) $$P=E_r{\text{Σ}}_f{\int }_0^a{\int }_{-H/2}^{H/2}A{J}_0\left(\frac{\mathrm{2,405}r}{a}\right)\times \mathrm{}cos\mathrm{}\left(\frac{\pi z}{H}\right)2\pi rdzdr,$$

Lembremos que o jacobiano em coordenadas cilíndricas é igual a $rdrd\theta dz$.

A integral em $z$ é facilmente calculada e resulta em $\frac{2H}{\pi }$.

Com isso, ficamos com

(32) $$P=4AH{E}_r{\text{Σ}}_f{\int }_0^a{J}_0\left(\frac{\mathrm{2,405}r}{a}\right)rdr.$$

Para calcular a integral em $r$ devemos lembrar que

$$J_0\left(B_{r}r\right)=\sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n}{{\left(n!\right)}^2}{\left(\frac{B_{r}r}{2}\right)}^{2n}.$$

Com isso, a integral em $r$ é dada por

$${\int }_0^a\sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n}{{\left(n!\right)}^2}{\left(\frac{B_{r}r}{2}\right)}^{2n}rdr.$$

Ou ainda

$${\displaystyle \sum }_{n=0}^{\infty }\frac{{(-1)}^n}{{(n!)}^2}{(\frac{B_r}{2})}^{2n}\underset{0}{\overset{a}{{\displaystyle \int }}}{r}^{2n+1}dr.$$

O que nos fornece

$$\sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n}{{\left(n!\right)}^2}{\left(\frac{B_r}{2}\right)}^{2n}\frac{a^{2n+2}}{2n+2}.$$

Podemos escrever a expressão anterior na seguinte forma

$$2\frac{a}{B_r}\sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n}{2n!\left(n+1\right)!}{\left(\frac{B_{r}a}{2}\right)}^{2n+1}.$$

Essa última expressão pode ser escrita em termos da função de Bessel de ordem 1, $J_1\left(B_{r}a\right)$, ou seja

$$\frac{a}{B_r}{J}_1\left(B_{r}a\right).$$

Lembrando que $B_r=\mathrm{2,405}/a$, a Eq. (32) fica dada por

$$P=\mathrm{1,6632}{a}^{2}AH{E}_r{\text{Σ}}_f{J}_1\left(\mathrm{2,405}\right).$$

Assim,

$$A=\frac{P}{\mathrm{1,6632}{a}^{2}H{E}_r{\text{Σ}}_f{J}_1\left(\mathrm{2,405}\right)}.$$

E assim, a solução final do problema é

(33) $$\varphi (r,z)=\frac{P}{\mathrm{1,6632}{a}^{2}H{E}_r{\text{Σ}}_f{J}_1\left(\mathrm{2,405}\right)}\times J_0\left(\frac{\mathrm{2,405}r}{a}\right)cos\left(\frac{\pi z}{H}\right).$$

Agora, vamos calcular a quantidade de nêutrons que escapam do reator a cada segundo. Para isso, precisamos calcular a densidade de corrente. O operador gradiente nem coordenadas cilíndricas é dado por

$\nabla =\widehat{r}\frac{\partial }{\partial r}+\widehat{\theta }\frac{1}{r}\frac{\partial }{\partial \theta }+\widehat{z}\frac{\partial }{\partial z}.$

Como para o caso considerado o fluxo de nêutrons não depende de $\theta$, temos que os componentes não-nulos da densidade de corrente de nêutrons são

(34) $$J_r=-D\frac{P}{\mathrm{1,6632}{a}^{2}H{E}_r{\text{Σ}}_f{J}_1\left(\mathrm{2,405}\right)}\times \left(\frac{\partial }{\partial r}{J}_0\left(\frac{\mathrm{2,405}r}{a}\right)\right)cos\left(\frac{\pi z}{H}\right).$$

e

(35) $$J_z=-D\frac{P}{\mathrm{1,6632}{a}^{2}H{E}_r{\text{Σ}}_f{J}_1\left(\mathrm{2,405}\right)}\times J_0\left(\frac{\mathrm{2,405}r}{a}\right)\left(\frac{\partial }{\partial z}cos\left(\frac{\pi z}{H}\right)\right).$$

Usando a relação $\frac{d{J}_0\left(x\right)}{dx}=-J_1\left(x\right)$, temos

(36) $$J_r=\frac{\mathrm{2,405}DP}{\mathrm{1,6632}{a}^{3}H{E}_r{\text{Σ}}_f{J}_1\left(\mathrm{2,405}\right)}\times J_1\left(\frac{\mathrm{2,405}r}{a}\right)cos\left(\frac{\pi z}{H}\right).$$

e ainda,

(37) $$J_z=\frac{\pi DP}{\mathrm{1,6632}{a}^2{H}^2{E}_r{\text{Σ}}_f{J}_1\left(\mathrm{2,405}\right)}\times J_0\left(\frac{\mathrm{2,405}r}{a}\right)\sin \left(\frac{\pi z}{H}\right).$$

Calculando o fluxo através das paredes do reator, ${\int }_A\overrightarrow{J}\cdot \widehat{n}dA$, obtemos

$${\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_{\raisebox{1ex}{$H$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}^{\raisebox{1ex}{$H$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{J}_r(r=a)ad\theta dz,}}$$

e

$$2{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^a{J}_z(z=H/2)rd\theta dr.}}$$

No cálculo do fluxo através da base e tampa do reator foi usada a relação $\int x{J}_0\left(x\right)dx=x{J}_1\left(x\right)$. Com isso, a quantidade de nêutrons que escapa por segundo através das paredes do reator é dada por

$$\frac{DP}{2a^2{H}^2{E}_r{\text{Σ}}_f}\left(\mathrm{11,57}{H}^2+{\pi }^2{a}^2\right)\text{nêutrons}/s.$$

Os resultados estabelecidos nesta seção serão comparados aos obtidos a partir das PINNs.

4. O Que é Uma Rede Neural Informada Por Física (PINN)?

Recentemente tivemos grandes avanços com as redes neurais para se fazer previões baseadas e dados de pré-treinamento. Aplicações em deep learning desde reconhecimento de imagens, identificação de padrões em grandes volumes de dados, análise de tendências com aplicações no mercado financeiro e, principalmente, o desenvolvimento de Large Language Models (LLMs) para facilitar a interface máquina-humano. Nos últimos anos, tem sido proposto o treinamento de redes neurais para a resolução de equações diferenciais parciais (EDP), através das redes neurais informadas pela física (PINNS em inglês), tanto para a resolução na forma direta quanto na forma indireta. A resolução de uma EDP na forma direta refere-se ao fato de, dadas as condiçoes de contorno, as soluções são previstas através da equação diferencial proposta com os coeficientes em suas formas funcionais conhecidos. O problema inverso refere-se a descobrir os coeficientes a partir de dados conhecidos.

Em ambos os casos, podemos destacar as vantagens em usarmos PINNs ao invés de outros métodos com um custo computacional maior como o método de elementos finitos (MEF):

1. Dispensam discretização espacial:

   • –

     É necessária uma malha computacional para o MEF, o que pode ser complexo em geometrias irregulares ou em altas dimensões.

   • –

     As PINNs aprendem a solução diretamente a partir das coordenadas espaciais (ou temporais), sem necessidade de discretização.

2. São mais indicadas para dados incompletos ou ruidosos:

   • –

     Os dados experimentais podem ser incorporados diretamente pelas PINNs, mesmo que esparsos ou ruidosos, ajustando a solução de forma híbrida (física + dados).

   • –

     O MEF depende fortemente de condições de contorno bem definidas.

3. Tratamento Natural de Problemas Inversos:

   • –

     O MEF resolve, na grande maioria dos casos, problemas diretos, enquanto as PINNs podem inferir parâmetros desconhecidos da equação (problemas inversos) de forma mais direta.

4. Comprovadamente, tem menor custo computacional em problemas de alta dimensionalidade:

   • –

     O MEF sofre com a “maldição da dimensionalidade” (ex.: problemas 3D+tempo exigem malhas extremamente finas).

   • –

     As PINNs aproximam a solução com uma rede neural, que pode ser mais eficiente em domínios complexos.

5. Suavidade da solução:

   • –

     As PINNs fornecem soluções infinitamente diferenciáveis (útil para calcular derivadas da solução).

   • –

     O MEF geralmente usa funções de bases contínuas por partes, limitando a suavidade.

6. Facilidade em problemas multifísicos e acoplados:

   • –

     As PINNs podem incorporar múltiplas leis físicas simultaneamente na função de perda.

   • –

     No MEF, modelos multifísicose exigem acoplamento complexo entre diferentes formulações.

Contudo, devemos destacar algumas desvantagens das PINNS. O treinamento pode ser demorado, pois redes neurais profundas exigem ajuste fino e podem ter convergência lenta. Uma outra desvantegem que merece destaque é a dificuldade em capturar soluções com descontinuidades, pois o MEF lida melhor com choques ou descontinuidades (ex.: escoamentos compressíveis). Porém, seu uso vem se ampliando e os modelos vêm sendo cada vez mais validados, constituindo uma nova área de desenvolvimento para possíveis aplicações Físicas e industriais.

4.1. Redes Neurais informadas pela Física(PINNs)

O termo “informadas pela Física” advém do fato de definirmos a taxa de perda (loss rate) a partir da EDP em questão. Em nosso caso, queremos resolver uma equação diferencial dada na forma:

(38) $$\{\begin{array}{cc}{\nabla }^{2}u\left(\mathbf{\text{r}}\right)+k^{2}u\left(\mathbf{\text{r}}\right)=0,\hfill & \\ u\left(\mathbf{\text{r}}\right)=0\mathrm{if}\mathbf{\text{r}}\in \partial \text{Ω},\hfill & \end{array}$$

que é uma equação diferencial parcial do tipo Helmholtz, no qual $B^2=\left(k_{\infty }-1\right)/L^2$ como já foi explicado anteriormente.

Agora, queremos aproximar a solução com uma rede neural. Para tanto, considere a rede neural diagramada pela Figura 2. A questão aqui é que queremos aproximar a solução tal que:

(39) $$u_{\psi }\left(\mathbf{\text{r}}\right)\approx u\left(\mathbf{\text{r}}\right),$$

Figura 2
Rede neural com 4 neurônios por camada escondida densa.

onde $u_{\psi }:{ℝ}^3\to ℝ$ denota uma função com parâmetros $\psi$ a ser treinada pela rede neural. Esta deve ser organizada em $C+1$ camadas de neurônios. A primeira camada ${𝒩}^{\left(0\right)}$ é a chamada camada de entrada que deve conter as variáveis, ou seja,

(40) $${𝒩}^{\left(0\right)}\left(\mathbf{\text{r}}\right)=\mathbf{\text{r}}=(x,y,z).$$

Cada camada $c$ subsequente será então parametrizada por um peso ${\mathbf{\text{P}}}^{\left(c\right)}\in {ℝ}^{s_{c-1}\times s_c}$ e por um bias${\mathbf{\text{b}}}^{\left(c\right)}\in {ℝ}^{s_c}$, onde $s_c$ é definido como tamanho da saida da camada $c$. As camadas $c\in [1,C-1]$ são chamadas de camadas escondidas densas cuja saida de cada uma é definida recursivamente como sendo:

(41) $${𝒩}^{\left(c\right)}\left(\mathbf{\text{r}}\right)=\xi \left({\mathbf{\text{P}}}^{\left(c\right)}{𝒩}^{\left(c-1\right)}\left(\mathbf{\text{r}}\right)+{\mathbf{\text{b}}}^{\left(c\right)}\right),$$

onde $\xi$ é chamada de função de ativação não-linear que especificaremos mais adiante. Ao final do processo, na camada de saída que terá apenas um valor será então dada por:

(42) $${𝒩}^{\left(C\right)}\left(\mathbf{\text{r}}\right)={\mathbf{\text{P}}}^{\left(C\right)}{𝒩}^{\left(C-1\right)}\left(\mathbf{\text{r}}\right)+{\mathbf{\text{b}}}^{\left(C\right)},$$

onde, finalmente, $u_{\psi }\left(\mathbf{\text{r}}\right)={𝒩}^{\left(C\right)}\left(\mathbf{\text{r}}\right)$.

Diagramaticamente, na Figura 2, temos a camada de entrada ${𝒩}^{\left(0\right)}$ representada por 3 neurônios para abrigar os valores de $x,y$ e $z$. Seguem então mais 3 camadas densas escondidas de 4 neurônios cada, finalizando na camada de saída ${𝒩}^{\left(C\right)}\left(\mathbf{\text{r}}\right)$. Neste caso, temos que $C=4$ e a chave é definirmos quem são as respectivas taxas de perda (loss rate), o que faremos logo abaixo. O número de parâmetros a serem treinados seguem a fórmula $p_e\ast E_1+E_1\ast E_2+E_2\ast E_3+E_3\ast p_s+E_1+E_2+E_3+p_s$. No caso da Figura 2 temos que $p_e=3$ (número de parâmetros de entrada), $E_1=E_2=E_3=4$ (número de neurônios por camada densa escondida) e $p_s=1$ (número de parâmetros de saída), totalizando 61 parâmetros a serem treinados pela rede neural. Este é apenas um exemplo esquemático mais simples para explicarmos o que faremos a seguir.

Na Figura 3 temos o dimensionamento real de parâmetros que usaremos em nossa PINN que serão $p_e=3$, $E_1=E_2=E_3=E_4=E_5=E_6=E_7=64$ e $p_s=1$, totalizando 25.281 parâmetros de treinamento da rede neural. A fim de capturarmos as particularidades da PDE em todos os pontos avaliados no cilindro nós utilizaremos 7 camadas densas de neurônios escondidos. Para tando, é necessário utlizarmos uma topologia de rede neural um pouco mais sofisticada para evitarmos um explosão de gradientes. Isso é alcançado criando-se “atalhos” através de add layers entre as camadas densas escondidas. Isso é feito automaticamente pelo tensorflow através da flag “PINN_ResNet” ao construirmos o modelo, como veremos a seguir no programa. O tensorflow mostra inicialmente a camada de entrada se conectando à primeira camada densa escondida como na Figura 2 e, em seguida, a concexão com as demais. Dessa forma, o tensorflow adicionará os add layers a cada duas camadas densas escondidas, criando os mencionados “atalhos” para os gradientes de aprendizado ao longo da rede. Vale ressaltar que esses add layers não contribuem com novos parâmetros de treinamento e são apenas vias de fluxo para os gradientes gerados entre cada duas camadas. Essa é uma característica fundamental distinta dessas redes, o que não ocorre com as redes sequenciais usuais como a da Figura 2.

Figura 3
Arquitetura utilizada para treinar a PINN do problema cilíndrico.

A técnica que empregaremos é a aproximação por aprendizado supervisionado usando dados de treinamento que, em nosso caso, serão os valores conhecidos da função em $\partial \text{Ω}$ e no ponto central. Seja então, $𝒟$ o conjunto com $N_{dados}$ desses dados, $(\mathbf{\text{r}},u)$, para treinamento:

(43) $$𝒟={\left\{({\mathbf{\text{r}}}_i^{dados},u_i^{dados})\right\}}_{i=1}^{N_{dados}}$$

Assim, $u_{\psi }\left(\mathbf{\text{r}}\right)$ será considerada uma boa aproximação se as previsões de $u$ estiverem próximas de $u_i^{dados}$ para cada $i\in [1,N_{dados}]$. Ou seja, queremos o conjunto de parâmetros otimizados tais que:

(44) $${\psi }^{\ast }=\underset{\psi }{\min }{ℒ}_{dados}\left(\psi \right)$$

com

(45) $$\begin{array}{l}{ℒ}_{dados}\left(\psi \right)=\frac{1}{N_{dados}}{\displaystyle \sum _{i=1}^{N_{dados}}{\left|u_{\psi }\left({\mathbf{\text{r}}}_i\right)-u_i^{dados}\right|}^2}\\ ={ℒ}_{centro}+{ℒ}_{faces}\end{array}$$

que é o erro quadrático médio que define a taxa de perda com relação aos dados de treinamento.

Ao resolvermos a Eq. (44) estaremos otimizando os parâmetros $\psi$ com relação aos dados de treinamento que, em nosso caso, serão os pontos $({\mathbf{\text{r}}}_i^{dados},u_i^{dados})$ em $\partial \text{Ω}$ e no ponto central, como mencionamos logo acima. Na verdade, para resolvermos esta equação de otimização utilizamos o algoritmo de gradiente descendente ou backpropagation. Esse algoritmo depende de uma diferenciação automática que é realizada sem problemas pelo Tensorflow sobre ${ℒ}_{dados}\left(\psi \right)$. Esta tarefa é realizada dentro de uma tolerância desejada através de:

(46) $${\psi }^{j+1}={\psi }^j-\alpha {▽}_{\psi }ℒ\left({\psi }^j\right)$$

com $ℒ={ℒ}_{dados}$ para a $j$-ésima iteração que chamamos de época. $\alpha$ é a taxa de aprendizado e deve ser otimizado para algumas rodadas iniciais. Se for muito pequeno demorará muito a convergir e se for grande não otimizará corretamente os parâmetros. Então, rodamos para um número pequeno de épocas e ajustamos o seu valor que geralmente fica em torno de 10^-3 a 10^-2, inicialmente. Vamos utilizar como otimizador o chamado adaptive moment estimation (ADAM) que ajusta a taxa de aprendizado dos parâmetros durante o treino, a partir de um valor inicial que estipulamos acima do valor final desejado. Essa técnica de resolução é chamada de Vanilla-PINN[¹⁸].

5. Soluções da Equação de Transporte de Nêutrons obtidas a partir das PINNs

Nesta seção, apresentamos a solução do reator cilíndrico via PINN; e para esse fim será usado o conteúdo discutido nas duas seções precedentes. Assim, queremos resolver a Eq. (18) por meio do treinamento de uma PINN que satisfaça certas condições de contorno para um cilindro de altura $H$ e raio $a$. Tal como fizemos anteriormente, definiremos um sistema de coordenadas e um fluxo de neutrons normalizados. Ao final dos cálculos teremos constantes multiplicativas que fazem o sistema reatornar às dimensões físicas reais, facilitando o treinamento das PINNs. Dessa forma, faremos as seguintes normalizações;

(47) $$\begin{array}{l}r=\sqrt{x^2+y^2}\to \frac{r}{a}\\ z\to \frac{z}{H}\\ u(r,z)=\frac{\varphi (r,z)}{\zeta }\end{array}$$

com

(48) $$\zeta =\frac{P}{\mathrm{1,632}{a}^{2}H{E}_R{\text{Σ}}_f{J}_1\left(\mathrm{2,405}\right)}.$$

Dessa forma, não é difícil ver que a Eq. (18) se torna uma equação similar à equação de Helmholtz Eq. (38)

(49) $$\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u}{\partial r}\right)+r_{aH}^2\frac{{\partial }^{2}u}{\partial z^2}+k^{2}u=0,$$

com $k=aB$ e $r_{aH}=a/H$. Temos agora uma equação diferencial normalizada e parametrizada pela razão entre o raio $a$ e a altura $H$ e o parâmetro $k^2={\left(aB\right)}^2=k_r^2+k_z^2$, com $k_r\approx \mathrm{2,405}$ e $k_z=\pi r_{aH}$.

No entanto, esta equação, tal como está, não é conveniente para o treino de uma PINN devido à singularidade em $r=0$. Isso traz instabilidades ao modelo que não permitem o treinamento da PINN de forma correta. Para contornarmos isso, vamos reescrever o operador diferencial desta equação em coordenadas cartesianas para obtermos

(50) $$\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+r_{aH}^2\frac{{\partial }^{2}u}{\partial z^2}+k^{2}u=0.$$

Assim, a Eq. (33), que representa a solução analítica será dada por:

(51) $$u(r,z)=u(x,y,z)=J_0\left(k_{r}r\right)\mathrm{}\cos \mathrm{}\left(\pi z\right)$$

com $r=\sqrt{x^2+y^2}$ e as condições de contorno dadas pela Eq. (38).

No listing 1, podemos ver as bibliotecas que serão necessárias para resolvermos a Eq. (38) utilizando as técnicas de aprendizado de máquina. A primeira biblioteca é a tensorflow que contém todas as ferramentas necessárias para o treinamento de modelos de redes neurais que serão utilizadas. A segunda é a biblioteca padrão para manipulações matemáticas e técnicas de vetorização chamada numpy que se integra perfeitamente com a biblioteca do tensorflow. A biblioteca scipy.special fornece as funções de Bessel utilizadas no cálculo. As demais bibliotecas matplotlib.pyplot e mpl_toolkits.mplot3d são bibliotecas de visualização dos resultados já bem estabelecidas do Python. Note que o parâmetro $k^2$, definido acima , entra como $ksq$ e $r_{aH}$ como $r$_aH e, além disso, definimos a precisão do modelo no tensorflow em “float64”. Temos também os pesos para a taxa de perda com relação à EDP em questão, um ponto central onde será colocada a fonte e às condições de contorno, respectivamente com wpde, wcentere wbc.

Listing 1
Bibliotecas a serem chamadas.

No listing 2 mostramos a função que gera os pontos do mesh para treinamento. Adotamos as coordenadas cilíndricas $(r,\theta ,z)$ e também mostramos a função que converte as coordenadas cilíndricas para cartesianas, pois utilizaremos com frequência essa transformação. Na Figura 4 mostramos os pontos do mesh gerados. Por padrão, a não ser que se queira aumentar o número de pontos do mesh, teremos: n_boundary=3000, n_domain=8000, o que significam 3000 pontos nas faces e 8000 pontos na região interna do cilindro. A lista boundary_values contém os valores da função $u$ nas fronteiras, ou seja, são os valores da condição de contorno sobre $u$ em todos os boundary_points. domain_points é uma lista com todos os pontos no interior do reator cilíndrico normalizado com $r\in \left[\mathrm{0,1}\right]$ e $z\in \left[-\mathrm{0.5,0.5}\right]$.

Listing 2
Gerador do mesh para treinamento.

Figura 4
mesh para o treinamento da PINN.

No listing 3 definimos a arquitetura de nossa rede neural e cada camada densa escondida é nomeada. Como vimos na Figura 3, utilizaremos 7 camadas com 64 neurônios cada. Essa parte é determinada por tentativa e erro. Devido ao ponto central, essa arquitetura se mostrou mais eficiente na convergência. Outras com menos hidden layers podem ser testadas aumentando-se o número de neurônios. Além disso, utilizaremos o ativador como sendo a função não linear “$\tanh$”, com o kernel inicializador glorot_normal, o que garante uma convergência mais rápida e estável, pois evita que os pesos iniciais impliquem uma explosão nos gradientes. Esse ativador é essencial para os cálculos uma vez que possui segunda derivada não nula, o que não acontece com ativadores lineares como “relu”. Isso se deve ao fato de estarmos lidando com operadores diferenciais de segunda ordem. Esta arquitetura está representada na Figura 3. Note que o InputLayer tem dimensão 3, embora estajamos interessados apenas nas coordenadas $r$ e $z$. Na verdade, a solução da Eq. (18) utilizou da simetria com relação à variável angular $\theta$. Dessa forma, teremos sempre o conjunto de variáveis $(r,\theta ,z)$ que serão importantes na hora de convertermos para coordenadas cartesianas $(x,y,z)$ quando necessário, como veremos adiante. Essa é uma boa etratégia que ficará mais clara nos próximos listings. Por fim, note a flag “PINN_ResNet” na definição do modelo ao fim desta rotina. Como comentamos anteriormente, ela será responsável pela inclusão de “atalhos” para os gradientes, evitando que haja uma explosão destes.

Listing 3
Arquitetura da rede neural.

No listing 4 definimos a a EDP a ser resolvida e a solução exata para comparação com o modelo. A implementação é feita através do decorador $@$, utilizando a função tensorflow.GradientTape. O modo “persistent” é utilizado para que numa mesma operação de memória sejam calculadas as derivadas segundas da função. Aqui utilizamos dois GradientTape para calcularmos as derivadas segundas. Isso foi necessário devido ao fato de termos que transformar para coordenadas cartesianas para evitarmos problemas com ponbtos mal definidos devido ao ângulo $\theta$, levando ao treinamento errático com pesos estagnados num ponto trivial. Como os textos explicativos no listing afirmam, o tensorflow calcula as derivadas em coordenadas cilíndricas internamente e reatorna os valores já convertidos. Dessa forma, o modelo é treinado em coordenadas cartesianas e a equação diferencial será resolvida integralmente em coordenadas cilíndricas.

Listing 4
Programação da EDP no tensorflow.

No listing 5 temos a definição da taxa de perda total. Introduzimos a fonte de nêutrons no ponto central no qual o valor de $u\left(r=0,\theta =0,z=0\right)=1$, que coincide em coordenadas cartesianas em $\left(x=0,y=0,z=0\right)$. A função compute_loss calcula a taxa total de perda do modelo e os reatorna com os respectivos pesos multiplicados. Essa é a função que será utilizada para otimizar o algoritmo e sua minimização é adotada como critério de convergência utilizando-se a Eq. (46). Um valor a partir de 10^-5 é considerado um bom critério para convergência e, como veremos, refletirá na reprodução dos resultados analíticos até uma certa acurácia.

Listing 5
Funções de perda e cálculos de erros e resíduos.

Diferentemente do treinamento das redes neurais tradicionais em que temos um conjunto de dados para treinamento e, após fazê-lo, utilizamos o modelo para fazer previsões, as PINNs utilizam a equação diferencial (ou equação física) para treinar ao mesmo tempo que faz previsões[⁹]. Isso só é possível porque a equação diferencial física fornece a lei que governa a EDP em todos os pontos do espaço em que queremos treinar a rede neural. Recentemente, houve uma explosão de trabalhos empregando PINNs para resolver EDPs em 2D na literatura[¹⁸]. Nota-se que ao empregar as EDPs normalizadas em módulo 1, o número de pontos de treinamento não precisa ser tão grande e a convergência ocorre relativamente de forma rápida. Em 3D, as coisas são diferentes. O fato de adicionarmos mais uma dimensão faz com que devamos adicionar muitos pontos para atingirmos certa acurácia. Contudo, realizar o treinamento ao mesmo tempo que as condições de contorno sejam fixadas nem sempre garante uma convergência computacional rápida e estável. Para tanto, devemos fazer um pré-treinamento apenas com as condições de contorno para que os parâmetros de treino reconheçam desde o início os valores da função explicitados pelas condições de contorno. Isso garante que a PINN convergirá em tempo hábil, como veremos a seguir. No listing 6 vemos a forma das funções de pré-treinamento e de treinamento de nossa rede neural. note a inclusão da função “tf.clip_by_global_norm” dentro do calculo do gradiente, o que garante a prevenção de “surtos” no gradiente. Além disso, veja que a perda total na otimização de pré-treino é definida como a perda nas faces com peso 100 mais a perda no ponto central com peso 10. Isso garante que os parâmetros iniciais para esses pontos já terão parâmetros pré-inicializados para que o modelo a ser treinado não comece tudo a partir do zero.

Listing 6
Funções de pré-treino e treinamento das PINNs.

Agora, vamos construir o modelo com os pontos nos locais certos, gerando a Figura 4. Para tanto, basta digitar os comandos do listing 7 logo abaixo. Temos como saída também um esquemático com as redes neurais geradas através da instância “tf.keras.utils.plot_model”, o que gera a Figura 3. Note o comando %matplotlib ipympl que é um comando “mágico”, típico do uso do jupyter notebook, para efeitos de visualizações interativas. Aconselhamos o uso do jupyter para as simulções pois facilitam na identificação de erros.

Listing 7
Para gerar o modelo inicial com os pontos distribuidos e gerar a figura 4.

Vamos então preparar a parte de visualização dos resultados. Temos no listing 8 a função “plot_results” na qual passamos o modelo treinado e o “loss_history” para visualização. As saidas estão em “resultados.png”, “resultados_2D.png” e “loss_rate.png”.

Listing 8
Visualização dos resultados após o treinamento.

Finalmente, no listing 9 fazemos o pré-treinamento e o treinamento da PINN para gerar os resultados e visualizá-los nos gráficos.

Listing 9
Pré-treinamento e treinamento da PINN e a geração das figuras.

Na Figura 5, mostramos a taxa de perda (“loss_rate”) do modelo treinado. Este é calculado na função “compute_loss” no listing 5 onde temos “loss_bc” para boundary conditions, “loss_center” para o ponto central e “loss_pde” para a equação diferencial parcial em questão. Essa função é chamada no listing 5 para ser otimizada pelo tensorflow.

Figura 5
Taxa de perda (“loss_rate”) para o ensaio do reator cilíndrico com 30000 Épocas.

Na Figura 6, mostramos os resultados convergidos da PINN treinada para o reator cilíndrico para o plano definido por $z=0$. Note a delineação das bordas no plano em $z=0$, bem como o formato da superfície em comparação com o resultado exato. Sem o pré-treinamento das condições de contorno, mesmo colocando-se um número muito grande de Épocas para treino, haveria deformações neste locais. A boa convergência pode ser visualizada através do painel que mostra a diferença de valores para um conjunto de pontos pré-determinados. Note que a diferença máxima ficou em torno de 1.18$\times {10}^{-3}$, o que indica uma convergência razoável. Para este ensaio utilizamos cerca de 30000 Épocas. Mostramos também nesta figura a projeção no plano $XY$ num mapa de superfície para vermos manchas quentes e frias de acordo com as diferenças com relação à solução exata. Vemos que o modelo parece mostrar boa concordância com os valores exatos próximo à região central e um erro um pouco maior nas regiões um pouco mais afastadas do centro, mas com boa concordância nas fronteiras.

Figura 6
Resultados após a convergência do modelo para o plano definido por $z=0$.

Quando normalizamos a equação diferencial do reator cilíndrico vimos que $k^2={\left(aB\right)}^2=k_r^2+k_z^2$, com $k_r\approx \mathrm{2,405}$ e $k_z=\pi r_{aH}$. Dessa forma, poderemos calcular o raio do reator cilíndrico para um $r_{aH}=0.5$, o que implica um cilindro com $H=2a$. Para tanto, se considerarmos um combustível típico de um reator comercial, a base de Urânio, a literatura nos fornece para uma usina de água leve – como um PWR/BWR [⁵], que é o caso típico de usinas comerciais – o valor de referência para o comprimento de difusão dos nêutrons térmicos na água é de cerca de $\mathrm{2,5}$ cm. Essa ordem de grandeza para o comprimento de difusão é corroborada em livros-texto de física de reatores nucleares que listam para água leve $L\approx \mathrm{2,54}$ cm[¹⁹, ²⁰]. Por outro lado, o fator de multiplicação infinito assume valores no intervalo $k_{\infty }=\mathrm{1,02}-\mathrm{1,05}$, com valor típico de 1,03[¹⁹, ²⁰]. Considerando esses valores e as Eqs.(21), (22) e (28), teremos para o reator cilíndrico em que $H=2a$ o raio $a=26$ cm e a altura $H=52$ cm. Há algumas limitações no modelo teórico aqui considerado se compararmos com as dimensões dos reatores reais, pois não incluímos em nossa análise ingredientes obrigatórios nos reatores reais tais como substâncias moderadoras (água pesada ou grafite); não levamos em consideração a heterogeneidade entre as células de combustível; no caso real, o fluxo costuma anular numa fronteira extrapolada, não exatamente na superfície geométrica como aparece na hipótese do nosso modelo. Essas nuances modificariam bastante os resultados obtidos para o tamanho do reator. Por esse motivo, essa situação idealizada é chamada de reator nu uniforme, mas possui elevada importância nas simulações reais para o entendimento do reator como um todo. A simplificação na análise se deve ao caráter pedagógico deste trabalho, ou seja, consideramos até o momento a equação de difusão de nêutrons aplicada a reatores nucleares na sua forma mais simples; assim, uma alternativa seria, além de reconsiderar as limitações mencionadas, utilizarmos, por exemplo, a equação de difusão de nêutrons multigrupos.

6. Considerações Finais e Perspectivas

Neste trabalho foram apresentados de forma pedagógica dois assuntos que despertam bastante interesse da comunidade acadêmica, quais sejam: reatores nucleares e inteligência artificial. Nesse processo, é esperado que o leitor, mesmo o menos familiarizado com tais conteúdos, compreenda a essência da física de nêutrons aplicada ao reator nuclear, bem como os príncipios básicos das redes neurais informadas por física. O material se destaca pela forma pedagógica que é apresentado, isto é, tanto as passagens matemáticas concernentes à equação de difusão de nêutrons, como a discussão acerca das PINNs foram construídas com bastante detalhe, com o objetivo de ser suficientemente inteligível. Nesse bojo, o problema do reator nuclear foi discutido analiticamente para o formato cilíndrico. o qual foi também solucionado com o uso da PINN, e os resultados obtidos por meio de ambos os métodos foram compatíveis. As limitações do modelo utilizado ficaram evidentes quando o tamanho estimado foi comparado com o tamanho real do reator; tendo tais diferenças origem nas hipóteses do aracbouço simplificado. Assim, como perspectiva, apontamos o aprofundamento, com o uso de PINNs, na equação de nêutrons multigrupo, assim seria possível tratarmos situações mais realísticas. Nessa perspectiva, o material apresentado apresenta potencial de ser utilizado por estudantes de graduação em física e engenharia que estejam interessados a iniciar estudos em física de reatores nucleares.

Agradecimento

Os autores agradecem ao Deputado Federal José Airton Félix Cirilo que através da emenda parlamentar N^o 23310003 tornou possível a rea-lização deste trabalho.

Disponibildade de dados

Todo o conjunto de dados que dá suporte ao presente texto foi publicado no corpo do artigo

Referências

• ^[1] G.I. Bell e S. Glasstone, Nuclear reactor theory (Van Nostrand Reinhold, New York, 1970).
• ^[2] J.J. Duderstadt e L.J. Hamilton, Nuclear reactor analysis (John Wiley & Sons, New York, 1976).
• ^[3] A.F. Henry, Nuclear-reactor analysis (MIT Press, Cambridge, 1975).
• ^[4] E.E. Lewis, Fundamentals of Nuclear Reactor Physics (Academic Press, San Dieg, 2008).
• ^[5] J.C. Lee, em: Wiley Encyclopedia of Electrical and Electronics Engineering (Wiley, Hoboken, 2016), 3 ed.
• ^[6] L. Vlahos, H. Isliker, Y. Kominis e K. Hizanidis, arXiv:0805.0419 (2008).
• ^[7] W.M. Stacey, Nuclear Reactor Physics (Wiley, New York, 2001).
• ^[8] J.R. Lamarsh e A.J. Baratta, Introduction to Nuclear Engineering (Prentice Hall, Upper Saddle River, 2001), 3 ed.
• ^[9] M. Raissi, P. Perdikaris e G.E. Karniadakis, J. Comp. Phys. 378, 686 (2019).
• ^[10] X. Jin, S. Cai, H. Li e G. Karniadakis, Journal of Computational Physics 426, 109951 (2021).
• ^[11] H. Jin, M. Mattheakis e P. Protopapas, em: 2022 International Joint Conference on Neural Networks– IJCNN (Piscataway, 2022).
• ^[12] G. Kissas, Y. Yang, E. Hwang, W.R. Witschey, J.A. Detre e P. Perdikaris, Computer Methods in Applied Mechanics and Engineering 358, 112623 (2020).
• ^[13] E. Schiassi, M. Florio, B.D. Ganapol, P. Picca e R. Furfaro, Annals of Nuclear Energy 167, 108833 (2022).
• ^[14] J. Wang, X. Peng, Z. Chen, B. Zhou, Y. Zhou e N. Zhou, Annals of Nuclear Energy 176, 109234 (2022).
• ^[15] Y. Yang, H. Gong, S. Zhang, Q. Yang, Z. Chen, Q. He e Q. Li, Annals of Nuclear Energy 183, 109656 (2023).
• ^[16] R.E. Marshak, Journal of Applied Physics 18, 185 (1947).
• ^[17] G.B. Arfken e H.J. Weber, Mathematical Methods for Physicists (Elsevier, Amsterdam/London, 2005), 6 ed.
• ^[18] H. Baty, arXiv:2403.00599 (2024).
• ^[19] J.A. Dejuren e H. Rosenwasser, Journal of Research of the National Bureau of Standards 51, 203 (1953).
• ^[20] M.M. Bretscher, Proceedings of the Indiana Academy of Science 72, 249 (1962).

Editado por

• Editor-Chefe:
  Marcello Ferreira https://orcid.org/0000-0003-4945-3169

Datas de Publicação

• Publicação nesta coleção
  14 Set 2026
• Data do Fascículo
  2026

Histórico

• Recebido
  08 Jun 2026
• Revisado
  31 Jul 2026
• Aceito
  02 Ago 2026