Apr-21-2023, 04:57 PM
Hey, so i`ve been working on this code for a while for my research but it just doesn't seem to work well, i think it is misscalculating in some step, probaly on the relative velocities, because it keeps exploding the value for some really high values, but i am not shure. Can someone help me finding the source of the problem here?
import numpy as np import matplotlib.pyplot as plt import pandas as pd import random as rd from matplotlib import animation from matplotlib.animation import PillowWriter # Constants kB = 1.38e-23 # Boltzmamn constant amu = 1.66e-27 # Atomic mass unity m_He = 4 * amu # He mass m_Ar = 38.9 * amu # Ar mass # Initial conditions n = 100000 # Total number os particles L = 10. # Length of the box T_He = 300 # initial temperature of He T_Ar = 300 #T_mistura = 500 #Temperature of the mixture conc_He = 0.90 # molar fraction of He in the mixture #T_Ar = (T_mistura - conc_He*T_He)/(1.-conc_He) #initial temperature of Ar T_mistura = conc_He*T_He + (1. - conc_He)*T_Ar n_He = int(n * conc_He) # Number of He particles n_Ar = n - n_He # Number of Ar particles sigma_He = np.sqrt((kB * T_He) / m_He) # standard deviation for He speeds sigma_Ar = np.sqrt((kB * T_Ar) / m_Ar) #standard deviation for Ar speeds T = 100 #Total time for the simulations dt = 0.00002 # time interval va = np.sqrt(2.*kB*T_mistura/(m_He + m_Ar)) #termal speed of the mixture def generate(): # Generating the initial speed and position of the particles x = L * np.random.rand(n, 3) # random positions v = np.zeros((n, 3),dtype='float64') # velocities array v[:n_He,:] = np.random.normal(0, sigma_He, size=(n_He, 3)) # He velocities v[n_He:,:] = np.random.normal(0, sigma_Ar, size=(n_Ar, 3)) # Ar velocities return v/va, x #va to let it dimentionless def plot(v): fig, ax = plt.subplots(1,2, figsize=(25,10)) a = ax[0] v_mag_He = np.linalg.norm(v[:n_He,:], axis=1) # Magnitude das velocidades do hélio a.hist(v_mag_He, bins=50, density=True, alpha=0.5, label='He') # Plot the histogram for the initial velocities of Argon v_mag_Ar = np.linalg.norm(v[n_He:,:], axis=1) # Magnitude das velocidades do argônio a.hist(v_mag_Ar, bins=50, density=True, alpha=0.5, label='Ar') # Configuring the graph a.set_title('Distribuição de velocidades') a.set_xlabel('Velocidade (m/s)') a.set_ylabel('Densidade de probabilidade') a.legend(fontsize=20) a.grid() a = ax[1] v_mag_He = v[:n_He,0] # Magnitude das velocidades do hélio a.hist(v_mag_He, bins=50, density=True, alpha=0.5, label='He') # Plot the histogram for the initial velocities of Argon v_mag_Ar = v[n_He:,0] # Magnitude das velocidades do argônio a.hist(v_mag_Ar, bins=50, density=True, alpha=0.5, label='Ar') # Configuring the graph a.set_title('Distribuição de velocidades') a.set_xlabel('Velocidade (m/s)') a.set_ylabel('Densidade de probabilidade') a.legend(fontsize=20) a.grid() plt.show() #Importing the matrix for viscosity muHA = pd.read_csv('tables/muHA.dat', sep="\s+",index_col=False).set_index(['0']) #Importing the matrix for cosx and the cross sections xiHe4He4 = pd.read_csv('tables/xiHe4He4.dat', sep="\s+", skiprows=4, header=None) xiHe4Ar = pd.read_csv('tables/xiHe4Ar.dat', sep="\s+", skiprows=4, header=None) xiArAr = pd.read_csv('tables/xiArAr.dat', sep="\s+", skiprows=4, header=None) def mover_particulas(v, x): # Move the particles in every Δt of the simulation #x = x + v*dt x,v = bordas(x,v) return x def temp_cinetica(vel): u = (1./n)*np.sum(v,axis=0) T_He = (m_He*va**2)/(3. * n_He * kB)* np.sum(v[:n_He]**2) T_Ar = (m_Ar*va**2)/(3. * n_Ar * kB)* np.sum(v[n_He:]**2) T = conc_He*T_He + (1.-conc_He)*T_Ar return T_He, T_Ar, T def bordas(pos, vel): """ Mantém as partículas dentro de uma caixa de tamanho L, fazendo com que elas refletam nas paredes da caixa quando colidem com elas. Args: pos (numpy array): Matriz Nx3 com as posições das N partículas. vel (numpy array): Matriz Nx3 com as velocidades das N partículas. L (float): Tamanho da caixa quadrada. Returns: numpy array: Matriz Nx3 com as novas posições das N partículas. numpy array: Matriz Nx3 com as novas velocidades das N partículas. """ # Verifica se alguma partícula está fora da caixa. new_pos = pos.copy() new_vel = vel.copy() # Verifica se as partículas estão dentro da caixa e corrige as posições e velocidades caso contrário for i in range(pos.shape[0]): for j in range(pos.shape[1]): if new_pos[i,j] < 0: new_pos[i,j] = abs(new_pos[i,j]) new_vel[i,j] = -new_vel[i,j] elif new_pos[i,j] > L: new_pos[i,j] = 2*L - new_pos[i,j] new_vel[i,j] = -new_vel[i,j] return new_pos, new_vel def colisao(v, tipo,sg): # type 1 = He-He # type 2 = Ar-Ar # type 3 = He-Ar p1 = 0 p2 = 0 while p1==p2: if tipo==1: #Selecting particles in the interval [0, n_He] p1 = int(np.floor((n_He)*rd.random())) p2 = int(np.floor((n_He)*rd.random())) #Particles masses m1 = m_He m2 = m_He #G constant for the type of collision G = 400. #Cross section matrix σ = xiHe4He4[100].to_numpy() #matrix for the values of cosX ξ = xiHe4He4.to_numpy() elif tipo == 2: #Selecting particles in the interval (n_He, n_Ar] p1 = int(n_He + np.floor(rd.random()*n_Ar)) p2 = int(n_He + np.floor(rd.random()*n_Ar)) #Particles masses m1 = m_Ar m2 = m_Ar #G constant for the type of collision G = 100. #Cross section matrix σ = xiArAr[100].to_numpy() ξ = xiArAr.to_numpy() else: #Selecting particles in the interval [0, n_He] and (n_He, n_Ar] p1 = int(np.floor((n_He)*rd.random())) p2 = int(n_He + np.floor(rd.random()*n_Ar)) #Particle masses m1 = m_He m2 = m_Ar #G constant for the type of collision G = 300. #Cross section matrix σ = xiHe4Ar[100].to_numpy() #matrix for the values of cosX ξ = xiHe4Ar.to_numpy() v1 = v[p1,:] #speed of the particle 1 v2 = v[p2,:] #speed of the particle 2 #Calculating relative speed g = v1 - v2 gr = np.linalg.norm(g) #norm of g #Post-collisional relative speed gl = np.zeros(3,dtype='float64') #Center of mass velocity G_cm = (m1*v1 + m2*v2)/(m1 + m2) #Calcuting the index of cross section array k = int(np.floor((np.log(1. + gr*va/G)/np.log(1.005)) + 0.5)) if k>899:k=899 vr = np.sqrt(g[1]*g[1] + g[2]*g[2]) if ((σ[k]*gr)/(sg) > rd.random()): #Accept-rejection condition print(σ[k], gr) n = int(np.floor(rd.random()*100)) #Monte Carlo parameters fot the post collisional velocities e = 2. * rd.random() * np.pi cosx = ξ[k][n] sinx = np.sqrt(1. - cosx*cosx) if vr>1e-15: #Components of the post collisional relative speed gl[0] = g[0]*cosx + vr*sinx*np.sin(e) gl[1] = g[1]*cosx + (gr*g[2]*np.cos(e) - g[0]*g[1]*np.sin(e))*sinx / vr gl[2] = g[2]*cosx - (gr*g[1]*np.cos(e) + g[0]*g[2]*np.sin(e))*sinx / vr else: gl[0] = g[0]*cosx gl[1] = g[0]*sinx*np.cos(e) gl[2]= g[0]*sinx*np.sin(e) #calculating the post collision velocities v1 = G_cm + 0.5*gl v2 = G_cm - 0.5*gl #setting the new velocities v[p1,:] = v1 v[p2,:] = v2 return v, σ[k]*gr else: return v, sg t=0. μ = muHA[f"{conc_He}"][np.floor(T_mistura)] σg11 = 15. σg12 = 15. σg22 = 15. v, x = generate() v_t = np.zeros((n, 3, T)) Temp = np.zeros((T,3)) for t in range(0,T,1): #plot(v) #x = mover_particulas(v,x) np.append(Temp[t,:],temp_cinetica(v)) N_col_HH = int((1./n)*n_He*(n_He - 1.)*μ*dt*σg11) N_col_AA = int((1./n)*n_Ar*(n_Ar - 1.)*μ*dt*σg22) N_col_HA = int((2./n)*n_He*n_Ar*μ*dt*σg12) print(f'N11 = {N_col_HH}, N12 = {N_col_HA}, N22 ={N_col_AA}') print(f'σ11 = {σg11}, σ12 = {σg12}, σ22 ={σg22}') for i in range(0,N_col_HH): v, sg11 = colisao(v, 1,σg11) if (sg11 > σg11): σg11 = sg11 for i in range(0,N_col_AA): v, sg22 = colisao(v, 2,σg22) if (sg22 > σg22): σg22 = sg22 for i in range(0,N_col_HA): v, sg12 = colisao(v, 3,σg12) if (sg12 > σg12): σg12 = sg12 print(t) plt.plot(np.arange(0,T,1), Temp[:,0],label='He') plt.plot(np.arange(0,T,1), Temp[:,1],label='Ar') plt.plot(np.arange(0,T,1), Temp[:,2], label='T') plt.legend() fig, ax = plt.subplots(1,1, figsize=(10,10)) def animate(i): ax.clear() v_mag_He = np.linalg.norm(v_t[:n_He,:,i], axis=1) # Magnitude das velocidades do hélio ax.hist(v_mag_He, bins=50, density=True, alpha=0.5, label='He') # Plot the histogram for the initial velocities of Argon v_mag_Ar = np.linalg.norm(v_t[n_He:,:,i], axis=1) # Magnitude das velocidades do argônio ax.hist(v_mag_Ar, bins=50, density=True, alpha=0.5,label='Ar') # Configuring the graph #ax.title('Distribuição de velocidades') ax.set_xlabel('Velocidade (m/s)',fontsize=20) ax.set_ylabel('Densidade de probabilidade',fontsize=20) #ax.set_xlim(0,10) #ax.set_ylim(0,2) ax.grid() ani = animation.FuncAnimation(fig, animate,frames=T,interval=10) ani.save('ani.gif',writer='pillow')