Untitled

import numpy as np
from scipy.linalg import lu_factor, lu_solve

def create_grid(R, I):
    dr = R / I
    edges = np.linspace(0, R, I + 1)
    centers = (edges[:-1] + edges[1:]) / 2
    return dr, centers, edges

def diffusion_setup(R, I, D, Sig_a, S, BC, geometry):
    dr, centers, edges = create_grid(R, I)
    A = np.zeros((I + 1, I + 1))
    b = np.zeros(I + 1)

    if geometry == 0:
        S_area = np.ones_like(edges)
        S_area[0] = 0
        V = np.ones_like(edges) * dr
    elif geometry == 1:
        S_area = 2 * np.pi * edges
        V = np.pi * (edges[1:I + 1]**2 - edges[0:I]**2)
    elif geometry == 2:
        S_area = 4 * np.pi * edges**2
        V = (4/3) * np.pi * (edges[1:I + 1]**3 - edges[0:I]**3)

    A[I, I] = (BC[0] / 2 + BC[1] / dr)
    A[I, I - 1] = (BC[0] / 2 - BC[1] / dr)

    Dplus = 0
    for i in range(I):
        r = centers[i]
        Dminus = Dplus
        Dplus = 2 * (D(r) * D(r + dr)) / (D(r) + D(r + dr))
        A[i, i] = (Dplus * S_area[i + 1]) / (dr * V[i]) + Sig_a(r)
        b[i] = S(r) * V[i]

        if i > 0:
            A[i, i - 1] = (-Dminus * S_area[i]) / (dr * V[i])
            A[i, i] += (Dminus * S_area[i]) / (dr * V[i])

        A[i, i + 1] = (-Dplus * S_area[i + 1]) / (dr * V[i])

    return centers, A, b

def solve_diffusion(R, I, D, Sig_a, S, BC, geometry):
    centers, A, b = diffusion_setup(R, I, D, Sig_a, S, BC, geometry)
    A_factored = lu_factor(A)
    phi = lu_solve(A_factored, b)
    return phi, centers

# Constants
radius_source = 2.3
S_source = 1e6
Pu_density = 19.74
Fe_density = 7.84
Pu_sigma_a = 2.11
Fe_sigma_a = 0.006
Pu_sigma_tr = 6.8
Fe_sigma_tr = 0.3703701
N_A = 6.022e23
M_Pu = 244
M_Fe = 56

# Number densities
N_Pu = (Pu_density * N_A) / M_Pu
N_Fe = (Fe_density * N_A) / M_Fe

# Macroscopic cross-sections
Sigma_a_Pu = N_Pu * Pu_sigma_a * 1e-24
Sigma_a_Fe = N_Fe * Fe_sigma_a * 1e-24
Sigma_tr_Pu = N_Pu * Pu_sigma_tr * 1e-24
Sigma_tr_Fe = N_Fe * Fe_sigma_tr * 1e-24

# Diffusion coefficients
D_Pu = 1 / (3 * Sigma_tr_Pu)
D_Fe = 1 / (3 * Sigma_tr_Fe)

# Functions for D, Sigma_a, and S
def D_func(r):
    return D_Pu if r < radius_source else D_Fe

def Sigma_a_func(r):
    return Sigma_a_Pu if r < radius_source else Sigma_a_Fe

def S_func(r):
    return S_source if r < radius_source else 0

# Boundary conditions
BC = [1, 0, 0]

# Geometry
geometry = 2

# Initial radius and cells
R = 10
I = 100

# Solve for flux
phi, centers = solve_diffusion(R, I, D_func, Sigma_a_func, S_func, BC, geometry)

# Calculate leakage
leakage = -D_func(R) * (phi[-1] - phi[-2]) / (centers[-1] - centers[-2])
print(f"Leakage: {leakage:.4e} neutrons/cm^2/s")

# Adjust shield thickness
while leakage > 1e4:
    R += 0.1
    phi, centers = solve_diffusion(R, I, D_func, Sigma_a_func, S_func, BC, geometry)
    leakage = -D_func(R) * (phi[-1] - phi[-2]) / (centers[-1] - centers[-2])
    print(f"Shield: {R:.2f} cm, Leakage: {leakage:.4e} neutrons/cm^2/s")

print(f"Shield thickness needed: {R - radius_source:.2f} cm")