Smoothed-Particle Hydrodynamics(SPH) in Python

What is Smoothed-particle hydrodynamics?

Smoothed-particle hydrodynamics(SPH)는 고체 역학 및 유체 흐름과 같은 유체 역학을 시뮬레이션 하는 데 사용되는 방법이다.

주로 참조한 논문은 「Smoothed Particle Hydrodynamics (SPH)를 이용한 유동 시뮬레이션 연구」(2015), 「Smoothed Particle Hydrodynamics」(The Gravitational Instability and its Role in the Evolution of Protostellar and Protoplanetary Discs, University of Leicester 2010) 이며, British Columbia 대학의 「Fluid Simulation」(2007) 교재의 도움을 많이 얻었다

SPH 기본 원리:

SPH 계산 과정:

SPH는 시뮬레이션 상황, 조건 등에 따라 다르게 구현되나, 일반적으로 아래의 단계를 거쳐 계산된다.

0. Smoothing Kernel Function

1. Calculate the density

2. Calculate the pressure

코드로 나타내면 아래와 같다.

def compute_density_and_pressure(particles, h):
    for i in range(len(particles)):
        particles[i].density = 0.0
        for j in range(len(particles)):
            if i != j:
                dist = np.sqrt((particles[i].x - particles[j].x)**2 + (particles[i].y - particles[j].y)**2)
                particles[i].density += particles[j].mass * poly6_kernel(dist, h)
        particles[i].pressure = gas_constant * (particles[i].density - rho_0)

3. Calculate the force

def compute_forces(particles, h):
    EPSILON = 1e-2
    for i in range(len(particles)):
        particles[i].fx = 0.0
        particles[i].fy = 0.0
        for j in range(len(particles)):
            if i != j:
                dist_x = particles[i].x - particles[j].x
                dist_y = particles[i].y - particles[j].y
                dist = np.sqrt(dist_x**2 + dist_y**2)

                # Softening Function to prevent division by zero
                softening_term = EPSILON * smoothing_length

                pressure_grad_x = -particles[i].mass * (particles[i].pressure + particles[j].pressure) / (
                    2 * (particles[j].density + softening_term)) * poly6_kernel(dist_x, h)
                pressure_grad_y = -particles[i].mass * (particles[i].pressure + particles[j].pressure) / (
                    2 * (particles[j].density + softening_term)) * poly6_kernel(dist_y, h)
        particles[i].fx += pressure_grad_x
        particles[i].fy += pressure_grad_y

4. Calculate the viscosity

def compute_vicosity_force(particles, h):
    EPSILON = 1e-2
    mu = dynamic_viscosity = 1.002
    for i in range(len(particles)):
        for j in range(len(particles)):
            if i != j:
                vel_x = particles[i].vx - particles[j].vx
                vel_y = particles[i].vy - particles[j].vy
                dist_x = particles[i].x - particles[j].x
                dist_y = particles[i].y - particles[j].y
                force_x = 0.0
                force_y = 0.0

                # Softening Function to prevent division by zero
                softening_term = EPSILON * smoothing_length

                force_x += particles[j].mass * (vel_x/(particles[j].density + softening_term)) * (poly6_kernel(dist_x, h)**2)
                force_y += particles[j].mass * (vel_y/(particles[j].density + softening_term)) * (poly6_kernel(dist_y, h)**2)

        particles[i].fx += mu * force_x
        particles[i].fy += mu * force_y

5. Add external forces and other equations

추가적으로 외부 힘들을 입자에 적용해준다.

여기서는 중력을 넣었다.

def update_gravity(particles, dt):
    for i in range(len(particles)):
        particles[i].vy += 9.8 * dt

입자의 이동을 더 질서 있게 만드는 XSPH correction을 추가로 적용했다. 공식은 다음과 같다.

def update_XSPH(particles, h):
    eplison = 1e-9
    EPSILON = 1e-9
    for i in range(len(particles)):
        for j in range(len(particles)):
            if i!=j:
                dist_x = particles[i].x - particles[j].x
                dist_y = particles[i].y - particles[j].y
                vel_x = particles[i].vx - particles[j].vx
                vel_y = particles[i].vy - particles[j].vy

                temp_correction_x = 0.0
                temp_correction_y = 0.0

                # Softening Function to prevent division by zero
                softening_term = EPSILON * smoothing_length

                temp_correction_x += (2*particles[j].mass * vel_x)/(particles[i].density+particles[j].density + softening_term)
                temp_correction_y += (2*particles[j].mass * vel_y)/(particles[i].density+particles[j].density + softening_term)
        temp_correction_x *= eplison * vel_x * poly6_kernel(dist_x, h)
        temp_correction_y *= eplison * vel_y * poly6_kernel(dist_y, h)

        particles[i].x += temp_correction_x
        particles[i].y += temp_correction_y

6. Update the acceleration and velocity

입자의 밀도와 압력으로부터 힘을 구했고, 해당 힘을 가지고 가속도, 속도, 위치를 업데이트 한다.

def update_particles(particles, dt):
    for i in range(len(particles)):
        particles[i].vx += particles[i].fx * dt / particles[i].mass
        particles[i].vy += particles[i].fy * dt / particles[i].mass
        particles[i].x += particles[i].vx * dt
        particles[i].y += particles[i].vy * dt

Particle Class

시뮬레이션을 위한 입자 클래스를 만들어 관리했다.

# Particle class
class Particle:
    def __init__(self, x, y, vx, vy, mass):
        self.x = x
        self.y = y
        self.vx = vx
        self.vy = vy
        self.fx = 0.0
        self.fy = 0.0
        self.mass = mass

전체 코드는 깃허브에서 확인할 수 있다.

https://github.com/bubbletok/IGB_SPH_Python

시뮬레이션 결과는 아래 동영상과 같다.

물 같은 유체를 표현하려고 했는데, 매우 ~~슬라임~~ 점도 높은 유체처럼 표현이 되었다.

추후에 spatial grid detection과 유체 방정식 등을 수정한 후 업데이트 할 예정이다.

References

[1] A Guided Journey into the Basics of the SPH Method

https://www.dive-solutions.de/blog/sph-basics

[2] Basics of the Smoothed Particle Hydrodynamics (SPH) Method

https://altair.com/newsroom/articles/Basics-of-the-Smoothed-Particle-Hydrodynamics-(SPH)-Method

[3] Smoothed-particle hydrodynamics

https://en.wikipedia.org/wiki/Smoothed-particle_hydrodynamics