particles_system.py

import numpy as np

k = 1.38062e-23  # Boltzmann constant


class ParticlesSystem:
    def __init__(self, radius: float, mass: float, volume: float, pressure: float, temperature: float,
                 init_type='norm'):
        """
        System of particles which interact like solid spheres
        :param radius: [m] radius of one particle
        :param mass: [kg] mass of one particle
        :param volume: [m^3] volume of system
        :param pressure: [pa] pressure in system
        :param temperature: [K] tempreture in system
        :param init_type: can take one of two values ['norm', 'uniform'] ('norm' by default);
                        'norm' - particles speeds will be generated by Maxwell-Boltzmann distribution
                        'uniform' - magnitude of velocities of each particle will be equal to mean square velocity
        """
        self.radius = radius
        self.mass = mass
        self.volume = volume
        self.pressure = pressure
        self.temperature = temperature
        self.concentration = pressure / (k * temperature)
        max_coordinate = self.volume ** (1 / 3) / 2
        step = self.concentration ** (-1 / 3)
        one_side = int(2 * max_coordinate // step)
        x = np.linspace(-max_coordinate, max_coordinate, one_side + 1)
        self.v_mean_square = np.sqrt(3 * k * temperature / mass)
        self.n_particles = int(np.round(one_side ** 3))
        self.p_coordinates = np.zeros((self.n_particles, 3), dtype=float)
        self.p_velocities = np.zeros((self.n_particles, 3), dtype=float)
        for i in range(one_side):
            for j in range(one_side):
                for g in range(one_side):
                    # initialization of particles positions
                    self.p_coordinates[one_side * one_side * i + one_side * j + g] = np.asarray([x[g], x[j], x[i]],
                                                                                                dtype=float)
                    if init_type == 'norm':
                        velocity = np.random.normal(loc=0, scale=np.sqrt(k * temperature / mass), size=3)
                    elif init_type == 'uniform':
                        direction = np.random.rand(3) - 0.5
                        direction = direction / np.linalg.norm(direction)
                        velocity = self.v_mean_square * direction
                    self.p_velocities[one_side * one_side * i + one_side * j + g] = velocity

    def iteration(self, dt):
        """
        Function iterates positions of particles in system and recalculates particle speeds after collisions
        :param dt: [sec] time of iteration
        :return:
        """
        if dt > 0.5 * self.radius / self.v_mean_square:
            print(f'Time of interaction ({0.5 * self.radius / self.v_mean_square}) is smaller than dt {dt}')
        # iteration of coordinates
        self.p_coordinates += self.p_velocities * dt
        # COLLISIONS of particle from p1 and corresponding particle from p2
        distances0 = self.p_coordinates[:, 0].reshape(-1, 1) - self.p_coordinates[:, 0].reshape(1, -1)
        distances1 = self.p_coordinates[:, 1].reshape(-1, 1) - self.p_coordinates[:, 1].reshape(1, -1)
        distances2 = self.p_coordinates[:, 2].reshape(-1, 1) - self.p_coordinates[:, 2].reshape(1, -1)
        distances = np.linalg.norm(np.stack((distances0, distances1, distances2)), axis=0)
        del distances0, distances1, distances2
        # finding colliding particles
        p1, p2 = np.where(distances <= 2 * self.radius)
        # deleting same particles from colliding pair of particles
        same_particles = np.where(np.equal(p1, p2) == True)[0]
        p1 = np.delete(p1, same_particles)
        p2 = np.delete(p2, same_particles)
        if p1.shape[0] != 0:
            parallel = (self.p_coordinates[p1] - self.p_coordinates[p2])
            parallel = parallel / np.linalg.norm(parallel, axis=1).reshape(-1, 1)
            v1_parallel = parallel * np.sum(self.p_velocities[p1] * parallel, axis=1).reshape(-1, 1)
            v1_norm = self.p_velocities[p1] - v1_parallel
            v2_parallel = parallel * np.sum(self.p_velocities[p2] * parallel, axis=1).reshape(-1, 1)
            v2_norm = self.p_velocities[p2] - v2_parallel
            v1_parallel_new = v2_parallel
            v2_parallel_new = v1_parallel
            self.p_velocities[p1] = v1_parallel_new + v1_norm
            self.p_velocities[p2] = v2_parallel_new + v2_norm
        # bounce from walls
        max_coordinate = self.volume ** (1 / 3) / 2 + self.radius
        particle, i = np.where(self.p_coordinates >= max_coordinate)
        self.p_velocities[particle, i] = -self.p_velocities[particle, i]
        particle, i = np.where(self.p_coordinates <= -max_coordinate)
        self.p_velocities[particle, i] = -self.p_velocities[particle, i]

    def get_velocities_magnitude(self):
        return np.linalg.norm(self.p_velocities, axis=1)

    def get_velocities_x(self):
        return self.p_velocities[:, 0]

    def get_velocities_y(self):
        return self.p_velocities[:, 1]

    def get_velocities_z(self):
        return self.p_velocities[:, 2]

# an example that can be used for initialization of system
# system = ParticlesSystem(radius=2e-10, mass=3.3e-26, volume=1e-23, pressure=10e5, temperature=300, init_type='norm')