Désintégration de l'uranium 235 : résolution numérique — Phénomènes de transport 3 : Diffusion de particules

L’uranium 235 n’a pas un noyau stable, celui-ci peut se fissionner en “captant” un neutron selon la réaction nucléaire

{\begin{matrix}  \end{matrix}}_{92}^{235} U

+ 1 neutron

\to

X

Y

𝜈

neutrons

où

X

Y

sont deux noyaux plus légers. La valeur moyenne de

𝜈

est

2,5

. Le nombre de réactions par unité de temps et de volume vaut

\frac{𝑛}{𝜏}

On se place en coordonnées sphériques, note

𝑛 (𝑡, 𝑟)

le nombre de neutrons par unité de volume et

\vec{𝑗} (𝑡, 𝑟)

le vecteur densité de courant de neutrons.

On prend pour condition aux limites

\forall 𝑡, 𝑛 (𝑡, 𝑟 = 𝑅) = 0

1/ En faisant un bilan de neutrons sur un volume mésoscopique, démontrer l’équation aux dérivées partielles vérifiée par

𝑛

Coup de pouce 1

Combien de neutrons sont captés durant

d 𝑡

dans le volume considéré ? Combien sont émis ?

Corrigé

On fait un bilan de neutrons sur une coquille sphérique de rayon intérieur $𝑟$ et d’épaisseur $d 𝑟$ . $d^{2} 𝑁 = 𝛿^{2} 𝑁$

d^{2} 𝑁 = \frac{𝜕 𝑛}{𝜕 𝑡} d 𝑉 d 𝑡

À chaque fois que la réaction se produit, $𝜈$ neutrons sont émis et un neutron est capté, donc $𝜈 - 1$ neutrons supplémentaires apparaissent. Le nombre de réactions dans le volume élémentaire durant $d 𝑡$ est $\frac{𝑛}{𝜏} d 𝑉 d 𝑡$ . Donc

\begin{aligned} 𝛿^{2} 𝑁 & = \frac{𝜈 - 1}{𝜏} 𝑛 d 𝑉 d 𝑡 + (𝑗 (𝑡, 𝑟) 4 𝜋 𝑟^{2} - 𝑗 (𝑡, 𝑟 + d 𝑟) 4 𝜋 {(𝑟 + d 𝑟)}^{2}) d 𝑡 \\ = \frac{𝜈 - 1}{𝜏} 𝑛 d 𝑉 d 𝑡 - 4 𝜋 \frac{𝜕}{𝜕 𝑟} (𝑟^{2} 𝑗 (𝑡, 𝑟)) d 𝑟 d 𝑡 \end{aligned}

Or $d 𝑉 = 4 𝜋 𝑟^{2} d 𝑟$ , donc en divisant par $d 𝑉 d 𝑡$ on trouve l’équation

\frac{𝜕 𝑛}{𝜕 𝑡} = \frac{𝜈 - 1}{𝜏} 𝑛 - \frac{1}{𝑟^{2}} \frac{𝜕}{𝜕 𝑟} (𝑟^{2} 𝑗 (𝑡, 𝑟))

La loi de Fick donne $\vec{𝑗} = - 𝐷 \vec{grad} 𝑛 = - 𝐷 \frac{𝜕 𝑛}{𝜕 𝑟} {\vec{𝑒}}_{𝑟}$ , avec $𝐷$ la diffusivité des neutrons dans l’uranium 235. En remplaçant dans l’équation, on trouve

\frac{𝜕 𝑛}{𝜕 𝑡} = \frac{𝐷}{𝑟^{2}} \frac{𝜕}{𝜕 𝑟} (𝑟^{2} \frac{𝜕 𝑛}{𝜕 𝑟}) + \frac{𝜈 - 1}{𝜏} 𝑛

2/ On pose

𝑦 (𝑡, 𝑟) = 𝑟 𝑛 (𝑡, 𝑟)

. Montrer que

𝑦

vérifie l’équation de diffusion

Coup de pouce 1

Remplacer

𝑛

par

\frac{𝑦}{𝑟}

dans l’équation obtenue précédemment.

Corrigé

On remplace $𝑛$ par $\frac{𝑦}{𝑟}$ dans l’équation précédente :

\begin{aligned} \frac{𝜕 \frac{𝑦}{𝑟}}{𝜕 𝑡} & = \frac{𝐷}{𝑟^{2}} \frac{𝜕}{𝜕 𝑟} (𝑟^{2} \frac{𝜕 \frac{𝑦}{𝑟}}{𝜕 𝑟}) + \frac{𝜈 - 1}{𝜏} \frac{𝑦}{𝑟} \\ = \frac{𝐷}{𝑟^{2}} \frac{𝜕}{𝜕 𝑟} (𝑟 \frac{𝜕 𝑦}{𝜕 𝑟} - 𝑦) + \frac{𝜈 - 1}{𝜏} \frac{𝑦}{𝑟} \\ = \frac{𝐷}{𝑟^{2}} (\frac{𝜕 𝑦}{𝜕 𝑟} + 𝑟 \frac{𝜕^{2} 𝑦}{𝜕 𝑟^{2}} - \frac{𝜕 𝑦}{𝜕 𝑟}) + \frac{𝜈 - 1}{𝜏} \frac{𝑦}{𝑟} \\ = \frac{𝐷}{𝑟} \frac{𝜕^{2} 𝑦}{𝜕 𝑟^{2}} + \frac{𝜈 - 1}{𝜏} \frac{𝑦}{𝑟} \end{aligned}

En multipliant par $𝑟$ , on trouve l’équation

\frac{𝜕 𝑦}{𝜕 𝑡} = 𝐷 \frac{𝜕^{2} 𝑦}{𝜕 𝑟^{2}} + \frac{𝜈 - 1}{𝜏} 𝑦

3/ Quelles sont les conditions aux limites vérifiées par

𝑦

𝑟 = 0

𝑟 = 𝑅

Coup de pouce 1

La densité particulaire doit rester finie en

𝑟 = 0

Corrigé

En $𝑟 = 0$ , la densité particulaire doit rester finie, donc $𝑦 (𝑡, 0) = 0$ .

En $𝑟 = 𝑅$ , on a la condition aux limites $𝑛 (𝑡, 𝑅) = 0$ , donc $𝑦 (𝑡, 𝑅) = 𝑅 𝑛 (𝑡, 𝑅) = 0$ .

Pour résoudre numériquement l’équation de diffusion, on discrétise l’espace avec un pas

d 𝑟

et le temps avec un pas

d 𝑡

. On note

𝑦_{𝑖, 𝑗} = 𝑦 (𝑖 d 𝑡, 𝑗 d 𝑟)

4/ Établir le schéma d’Euler explicite pour résoudre numériquement l’équation de diffusion vérifiée par

𝑦

Coup de pouce 1

Utiliser la relation de Taylor à l’ordre 2 pour approximer

𝑦_{𝑖, 𝑗 + 1}

𝑦_{𝑖, 𝑗 - 1}

. De même, approximer

𝑦_{𝑖 + 1, 𝑗}

à l’ordre 1.

Corrigé

$𝑦_{𝑖, 𝑗 + 1} = 𝑦 (𝑖 d 𝑡, 𝑗 d 𝑟 + d 𝑟) \approx 𝑦 (𝑖 d 𝑡, 𝑗 d 𝑟) + \frac{𝜕 𝑦}{𝜕 𝑟} d 𝑟 + \frac{1}{2} \frac{𝜕^{2} 𝑦}{𝜕 𝑟^{2}} {d 𝑟}^{2}$

$𝑦_{𝑖, 𝑗 - 1} = 𝑦 (𝑖 d 𝑡, 𝑗 d 𝑟 - d 𝑟) \approx 𝑦 (𝑖 d 𝑡, 𝑗 d 𝑟) - \frac{𝜕 𝑦}{𝜕 𝑟} d 𝑟 + \frac{1}{2} \frac{𝜕^{2} 𝑦}{𝜕 𝑟^{2}} {d 𝑟}^{2}$

En additionnant, on trouve $𝑦_{𝑖, 𝑗 + 1} - 2 𝑦_{𝑖, 𝑗} + 𝑦_{𝑖, 𝑗 - 1} \approx \frac{𝜕^{2} 𝑦}{𝜕 𝑟^{2}} {d 𝑟}^{2}$

De plus, $𝑦_{𝑖 + 1, 𝑗} = 𝑦 (𝑖 d 𝑡 + d 𝑡, 𝑗 d 𝑟) \approx 𝑦 (𝑖 d 𝑡, 𝑗 d 𝑟) + \frac{𝜕 𝑦}{𝜕 𝑡} d 𝑡 = 𝑦_{𝑖, 𝑗} + \frac{𝜕 𝑦}{𝜕 𝑡} d 𝑡$ , d’où $\frac{𝜕 𝑦}{𝜕 𝑡} \approx \frac{𝑦_{𝑖 + 1, 𝑗} - 𝑦_{𝑖, 𝑗}}{d 𝑡}$ .

L’équation de diffusion s’écrit donc

\frac{𝑦_{𝑖 + 1, 𝑗} - 𝑦_{𝑖, 𝑗}}{d 𝑡} = 𝐷 \frac{𝑦_{𝑖, 𝑗 + 1} - 2 𝑦_{𝑖, 𝑗} + 𝑦_{𝑖, 𝑗 - 1}}{{d 𝑟}^{2}} + \frac{𝜈 - 1}{𝜏} 𝑦_{𝑖, 𝑗}

On peut isoler $𝑦_{𝑖 + 1, 𝑗}$ :

𝑦_{𝑖 + 1, 𝑗} = 𝑦_{𝑖, 𝑗} + 𝐷 \frac{d 𝑡}{{d 𝑟}^{2}} (𝑦_{𝑖, 𝑗 + 1} - 2 𝑦_{𝑖, 𝑗} + 𝑦_{𝑖, 𝑗 - 1}) + \frac{𝜈 - 1}{𝜏} d 𝑡 𝑦_{𝑖, 𝑗}

5/ Compléter le code Python ci-dessous pour simuler la désintégration de l’uranium 235 dans la boule. On prendra comme condition initiale une densité de neutrons uniforme dans la boule égale à

1

⋅1012⁠

⁠

m^{- 3}

(sauf aux conditions aux limites où elle est nulle).

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

import numpy as np

R = 0.2 # rayon de la boule en m
nu = 2.5 # nombre moyen de neutrons émis par fission
tau = 5.4e-9 # temps moyen entre deux fissions en s
D = 2e5 # diffusivité des neutrons dans l'uranium 235 en m2/s

Nx = 50 # nombre d'échantillons spatiaux
Nt = 1000 # nombre d'échantillons temporels

dx = ... # pas spatial
dt = 0.5 * dx**2 / (2*D) # pas temporel (condition de stabilité)

t = ... # array contenant tous les instants
r = ... # array contenant toutes les positions

y = ... # initialisation de la matrice avec des zéros
y[0,1:-1] = ... # condition initiale : densité uniforme 

for i in range(0, Nt-1):
    y[i+1, ...] = 0 # condition aux limites en r=0
    y[i+1, ...] = 0 # condition aux limites en r=R
    for j in range(1, Nx-1):
        # schéma d'Euler explicite
        y[i+1,j] = ...

Coup de pouce 1

On peut utiliser les fonctions np.linspace et np.zeros de la bibliothèque numpy.

Corrigé

1
2
3
4
5
6
7
8
9
10
11
12
13
14

dx = R / (Nx - 1) # pas spatial

t = np.arange(Nt) * dt # array contenant tous les instants
r = np.linspace(0, R, Nx) # array contenant toutes les positions

y = np.zeros( (Nt,Nx) ) # initialisation de la matrice avec des zéros
y[0,1:-1] = 1e12 * r[1:-1] # condition initiale : densité uniforme 

for i in range(0, Nt-1):
    y[i+1, 0] = 0 # condition aux limites en r=0
    y[i+1, -1] = 0 # condition aux limites en r=R
    for j in range(1, Nx-1):
        # schéma d'Euler explicite
        y[i+1,j] = y[i,j] + D*dt/dx**2 * (y[i,j+1] - 2*y[i,j] + y[i,j-1]) + (nu-1)/tau * y[i,j] * dt

6/ Tracer la densité de neutrons en fonction du rayon aux instants

𝑡 = 0

𝑡 = \frac{𝑡_{max}}{4}

𝑡 = \frac{𝑡_{max}}{2}

𝑡 = 3 \frac{𝑡_{max}}{4}

𝑡 = 𝑡_{max}

où

𝑡_{max}

est le temps final de la simulation.

Coup de pouce 1

Utiliser la bibliothèque matplotlib pour tracer les graphiques demandés.

Coup de pouce 2

Penser à calculer la densité de neutrons

𝑛

à partir de

𝑦

en utilisant la relation

𝑛 = \frac{𝑦}{𝑟}

Corrigé

On peut utiliser le code suivant pour tracer la densité de neutrons en fonction du rayon à différents instants :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

n = np.zeros_like(y)
n[:, 1:] = y[:, 1:] / r[1:] # calcul de la densité de neutrons (r[0] = 0)
n[:, 0] = n[:, 1] # en r=0, n est finie : on prolonge par continuité

import matplotlib.pyplot as plt
plt.figure()
plt.plot(r, n[0,:], label=f't={t[0]:.2e} s')
plt.plot(r, n[Nt//4,:], label=f't={t[Nt//4]:.2e} s')
plt.plot(r, n[Nt//2,:], label=f't={t[Nt//2]:.2e} s')
plt.plot(r, n[3*Nt//4,:], label=f't={t[3*Nt//4]:.2e} s')
plt.plot(r, n[-1,:], label=f't={t[-1]:.2e} s')
plt.xlabel('Rayon r (m)')
plt.ylabel('Densité de neutrons (m$^{-3}$)')
plt.title("Évolution de la densité de neutrons dans la boule d'uranium 235")
plt.legend()
plt.show()

7/ Tracer l’évolution temporelle de la densité de neutrons en

𝑟 = \frac{𝑅}{2}

Corrigé

1
2
3
4
5
6

plt.figure()
plt.plot(t, n[:, Nx//2])
plt.xlabel('Temps t (s)')
plt.ylabel('Densité de neutrons (m$^{-3}$)')
plt.title("Évolution temporelle de la densité de neutrons en r=R/2")
plt.show()

8/ Pour de petites valeurs de

𝑅

, la densité de neutrons tend vers

0

avec le temps. Pour de grandes valeurs de

𝑅

, la densité de neutrons croît exponentiellement avec le temps. Déterminer la valeur critique de

𝑅

séparant ces deux comportements. On pourra procéder par essais successifs et on la déterminera à

0,5

⁠

cm

près.

Coup de pouce 1

Le profil initial n’est pas le mode fondamental : il commence par se déformer, ce qui masque la croissance ou la décroissance exponentielle. Il faut donc simuler assez longtemps.

Corrigé

Avec $𝑁_{𝑡} = 1000$ , la simulation ne dure que quelques nanosecondes : le profil initial $𝑦 = 𝑛_{0} 𝑟$ n’est pas le mode fondamental $sin (𝜋 𝑟 / 𝑅)$ , il se déforme d’abord, et cette relaxation masque complètement la tendance exponentielle près du rayon critique. Il faut allonger la simulation, par exemple Nt = 20000.

Pour $𝑅 = 0,08 m$ , la densité de neutrons tend vers $0$ avec le temps.

Pour $𝑅 = 0,09 m$ , la densité de neutrons croît exponentiellement avec le temps.

Le rayon critique est situé entre ces deux valeurs : $𝑅_{𝑐} = (8,5 \pm 0, 5) cm$ , en accord avec la valeur analytique $𝜋 \sqrt{𝐷 𝜏 / (𝜈 - 1)} = 8,4 cm$ .

🖥️ Désintégration de l’uranium 235 : résolution numérique ★