Untitled
unknown
plain_text
6 months ago
5.5 kB
5
Indexable
import numpy as np import matplotlib.pyplot as plt import pandas as pd pd.set_option('display.max_rows', None) pd.set_option('display.max_columns', None) class NeutronFluxSimulation: def __init__(self, E_a, D, R_1, s0, n): self.E_a = E_a # cm^-1 self.D = D # cm self.R_1 = R_1 # cm self.s0 = s0 # Source strength self.n = n # Number of nodes self.r = np.linspace(0, self.R_1, self.n + 1) self.dr = self.r[1] - self.r[0] self.A = np.zeros((self.n + 1, self.n + 1)) self.b = np.zeros(self.n + 1) def setup_matrices(self): # Left boundary r_half = (self.r[0] + self.dr / 2) r_n_min = self.r[self.n] - (self.dr / 2) self.A[0, 0] = ((self.D * r_half) / self.dr) + (self.E_a * 0.5 * np.square(r_half)) self.A[0, 1] = (-self.D * r_half) / self.dr self.b[0] = self.s0 * 0.5 * np.square(r_half) # Right boundary self.A[self.n, self.n] = (self.r[self.n] / 2) + ((self.D * r_n_min) / self.dr) + ( self.E_a * 0.5 * (np.square(self.r[self.n]) - np.square(r_n_min))) self.A[self.n, self.n - 1] = (-self.D * r_n_min) / self.dr self.b[self.n] = self.s0 * 0.5 * (np.square(self.r[self.n]) - np.square(r_n_min)) # Interior nodes for i in range(1, self.n): # from 1 to n-1 r_i_min = self.r[i] - (self.dr / 2) r_i_plus = self.r[i] + (self.dr / 2) self.A[i, i] = ((self.D * (r_i_plus + r_i_min) / self.dr)) + ( self.E_a * 0.5 * (np.square(r_i_plus) - np.square(r_i_min))) self.A[i, i - 1] = -self.D * r_i_min / self.dr self.A[i, i + 1] = -self.D * r_i_plus / self.dr self.b[i] = self.s0 * 0.5 * (np.square(r_i_plus) - np.square(r_i_min)) def gauss_seidel(self, tol=1e-6, max_iterations=100000): n = len(self.b) x = np.zeros_like(self.b) for iteration in range(max_iterations): x_old = np.copy(x) for i in range(n): sum1 = np.dot(self.A[i, :i], x[:i]) sum2 = np.dot(self.A[i, i + 1:], x_old[i + 1:]) x[i] = (self.b[i] - sum1 - sum2) / self.A[i, i] # Check for convergence if np.linalg.norm(abs(x - x_old), ord=np.inf) < tol: print(f'Gauss-Seidel: Converged after {iteration + 1} iterations') return x print('Max iterations reached without convergence') return x def jacobi(self, tol=1e-6, max_iterations=100000): n = len(self.b) x = np.zeros_like(self.b) x_new = np.zeros_like(self.b) for iteration in range(max_iterations): for i in range(n): sum1 = np.dot(self.A[i, :i], x[:i]) sum2 = np.dot(self.A[i, i + 1:], x[i + 1:]) x_new[i] = (self.b[i] - sum1 - sum2) / self.A[i, i] # Check for convergence if np.linalg.norm(x_new - x, ord=np.inf) < tol: print(f'Jacobi: Converged after {iteration + 1} iterations') return x_new x = np.copy(x_new) print('Max iterations reached without convergence') return x_new def calculate_rates(self): fi_gauss = self.gauss_seidel() fi_jacobi = self.jacobi() EaM = np.zeros(self.n + 1) EaM[0] = self.E_a * np.pi * (self.r[0] + self.dr / 2) ** 2 * fi_gauss[0] EaM[self.n] = self.E_a * np.pi * fi_gauss[self.n] * (self.r[self.n] ** 2 - (self.r[self.n] - self.dr / 2) ** 2) for i in range(1, self.n): EaM[i] = self.E_a * np.pi * ((self.r[i] + self.dr / 2) ** 2 - (self.r[i] - self.dr / 2) ** 2) * fi_gauss[i] absorption_rate = np.sum(EaM) # Calculate the leakage rate at the boundary leakage_rate = (fi_gauss[self.n] / 2) * 2 * np.pi * self.R_1 source_rate = self.s0 * np.pi * self.R_1 ** 2 ab_leak = absorption_rate + leakage_rate data = {"Leakage Rate": [leakage_rate], "Absorption Rate": [absorption_rate], "Source Rate": [source_rate], "Absorption Rate + Leakage Rate": [ab_leak] } index_labels = ["Rate Values"] df = pd.DataFrame(data, index=index_labels) print(df) return fi_gauss, fi_jacobi def plot_results(self, fi_gauss, fi_jacobi): plt.figure(figsize=(10, 6)) plt.plot(self.r[:(self.n + 1)], fi_gauss[:(self.n + 1)], label='Neutron Flux (Gauss-Seidel)', color='blue') plt.xlabel('x (cm)') plt.ylabel('Flux (fi)', labelpad=15) plt.title('Neutron Flux Distribution using Gauss-Seidel Method') plt.grid(True) plt.legend() plt.tight_layout() plt.show() plt.figure(figsize=(10, 6)) plt.plot(self.r[:(self.n + 1)], fi_jacobi[:(self.n + 1)], label='Neutron Flux (Jacobi)', color='black') plt.xlabel('x (cm)') plt.ylabel('Flux (fi)', labelpad=15) plt.title('Neutron Flux Distribution using Jacobi Method') plt.grid(True) plt.legend() plt.tight_layout() plt.show() simulation = NeutronFluxSimulation(E_a=0.1, D=1.2, R_1=40, s0=10 ** 4, n=50) simulation.setup_matrices() fi_gauss, fi_jacobi = simulation.calculate_rates() simulation.plot_results(fi_gauss, fi_jacobi)
Editor is loading...
Leave a Comment