비압축성 유체의 효율적인 시뮬레이션은 여전히 어렵고 해결되지 않은 문제입니다.

기존 연구들은 종종 비압축성, 안정성, 그리고 비용 사이에서 다양한 절충을 합니다.

하지만, 이 세 가지를 모두 얻는 경우는 드뭅니다.

본 논문에서는 비압축성에 접근하기 위해 특별히 공식화된 변분 에너지를 최적화하고자 2차 암시적 하강 기법을 사용하는 새로운 비압축성 평활 입자 유체 동역학(SPH) 기법을 소개합니다.

저희는 제안하는 방법이 최소한의 계산 예산 비용

더 나아가, 저희 방법은 극한의 타임스텝 하에서도 무조건적으로 안정적이어서, 상호작용 애플리케이션에 적합하다는 것을 보여줍니다.

라그랑주 입자(Lagrangian particles) 기반의 smoothed particle hydrodynamics (SPH) Müller/2003/Particle-based fluid simulation는 컴퓨터 그래픽스에서 유체 시뮬레이션을 위한 대중적인 기법이었습니다.

이 기법은 종종 액체 상태의 물과 같은 비압축성 유체를 시뮬레이션하는 데 사용되지만, SPH로 비압축성을 달성하는 것은 근본적인 과제입니다.

이는 SPH 공식에서 비압축성이 유체 밀도를 유지하는 강한 제약 조건에 해당하기 때문입니다.

결과적으로, 명시적 속도 적분(explicit velocity integration) Müller/2003/Particle-based fluid simulation은 안정성을 위해 극도로 작은 타임스텝을 요구하며, 비선형성이 높은 제약 조건을 사용하여 비압축성을 달성하기 위해서는 훨씬 더 작은 타임스텝이 필요합니다 Becker/2007/Weakly compressible SPH.

따라서 대부분의 최근 연구들은 입자별 제약 조건 강성을 적응적으로 조절하는 준-암시적 속도 적분(semi-implicit velocity integrations)의 변형들을 채택해왔습니다 Bender/2015/Divergence-free SPH, Ihmsen/2014/Implicit Incompressible SPH.

불행히도, 이러한 접근법들은 복잡한 시뮬레이션에서 안정성 문제나 과도한 압축으로 여전히 어려움을 겪으며, 비교적 작은 타임스텝을 요구합니다.

이는 병렬 컴퓨팅을 통한 계산 성능이라는 SPH 기반 방법들의 가장 강력한 장점 중 하나를 저해합니다.

한 가지 주목할 만한 예외는 position-based fluids (PBF) 방법 Macklin/2013/Position based fluids으로, SPH를 위치 기반 동역학 Macklin/2016/XPBD, Müller/2007/Position based dynamics과 유사한 제약 조건 공식으로 변환하여 더 큰 타임스텝에서도 뛰어난 안정성을 제공합니다.

하지만, PBF에서 사용되는 1차 근사(first-order approximations)는 애니메이션의 품질을 저하시키고, 움직임을 상당히 감쇠시키며 물리적으로 타당한 결과에서 벗어나게 합니다.

본 논문에서는 SPH를 위한 또 다른 위치 기반 공식을 소개하며, 암시적 오일러의 변형 형태(variational form of implicit Euler)에서 시작하여 SPH에 대한 완전한 암시적 적분(fully implicit integration)을 제공합니다.

따라서 저희는 이 방법을 implicit position-based fluids (IPBF)라고 부르지만, 저희의 적분 방식은 PBF와 상당히 다릅니다.

저희는 vertex block descent (VBD) Chen/2024/Vertex Block Descent와 유사한 2차 하강 공식(second-order descent formulation)을 사용한 완화된 야코비 반복(relaxed Jacobi iterations)을 사용하여, 과도한 감쇠를 피하면서 무조건적인 안정성을 달성하고 더 큰 타임스텝을 허용합니다.

또한, 저희는 임의로 강한 밀도 제약 조건을 사용하여 비압축성을 모델링할 수 있으며, 더 큰 타임스텝에서 개선된 비압축성을 달성합니다.

강성 시스템(Stiff systems)은 특히 수렴을 방해하는 제한된 반복 횟수로 시뮬레이션될 때 에너지를 주입하는 것으로 알려져 있습니다.

IPBF도 예외는 아닙니다.

비록 무한한 강성에서도 에너지 주입은 비교적 미미하고 불안정성을 초래하지는 않지만, 이는 유체가 정지 상태에 도달하는 것을 막아 일부 응용 프로그램에서는 바람직하지 않을 수 있습니다.

따라서 저희는 덜 강한 제약 조건의 해와 비교하여 운동 에너지를 추출하는 간단한 감쇠 공식(simple damping formulation)도 소개합니다.

이 선택적 감쇠는 움직임이 잦아들 때까지 애니메이션에 미치는 영향을 최소화하도록 설계되었습니다.

저희의 실험은 동일한 계산 예산 하에서 저희 방법이 속도 기반 방법들보다 우수한 수치적 안정성을 제공하고, PBF를 포함한 이전 SPH 기반 방법들에 비해 개선된 비압축성 거동을 보인다는 것을 보여줍니다.

특히, 저희는 유체 움직임을 감쇠시키거나 부자연스러운 애니메이션을 생성하지 않으면서 이러한 성과를 달성합니다.

SPH는 원래 항성 간 유동을 시뮬레이션하기 위해 고안된 방법으로(Gingold and Monaghan/1977/SPH; Lucy/1977/SPH), 유체를 질량, 운동량 및 기타 필드 양을 저장하는 입자로 이산화합니다.

SPH는 Müller et al./2003/Particle-based fluid에 의해 컴퓨터 그래픽스 응용 분야에 도입되었으며, 이후 Koschier et al./2022/A Survey on SPH Methods에 기술된 바와 같이 SPH 기반 솔버에서 사용할 수 있는 추정 기법과 힘 계산을 개선하기 위한 방대한 연구로 이어졌습니다. 특히 압력 계산은 유체의 비압축성을 유지하는 것이 현실적이고 역동적인 시뮬레이션으로 이어지기 때문에 많은 연구의 주제가 되어 왔습니다.

표준 SPH(Müller et al./2003/Particle-based fluid)는 비압축성 유체를 위한 상태 방정식 압력 해법을 도입했는데, 여기서 각 입자는 정지 상태로부터의 밀도 편차를 보상하기 위해 페널티 힘을 계산합니다. 이후 비선형 공식으로 개선되었지만(Becker and Teschner/2007/Weakly compressible SPH), 이러한 상태 방정식 솔버는 현실적인 유체 흐름을 달성하기 위해 높은 강성 상수를 요구하며, 이로 인해 솔버가 수치적 안정성을 유지하기 위해 낮은 타임스텝을 사용하도록 강제합니다.

더 높은 타임스텝에서 안정성을 유지하기 위한 개선은 비압력 힘으로 인한 밀도 변동에 기반하여 페널티 힘을 계산하기 위해 예측-수정 기법을 사용하는 반복 솔버에 의해 이루어졌습니다(He et al./2012/Local Poisson SPH; Solenthaler and Pajarola/2009/Predictive-corrective incompressible SPH). Implicit Incompressible SPH (IISPH)와 같은 일부 솔버는 사용자가 정의한 강성 대신, 업데이트된 압력 힘을 계산할 때 입자별 강성 매개변수를 계산합니다(Ihmsen et al./2014/Implicit Incompressible SPH). 더 높은 타임스텝에서 안정성을 향상시켰지만, 이러한 방법들은 최소한의 반복 횟수와 큰 잔여 오차 누적이 주어졌을 때 비압축성을 유지하는 데 실패합니다.

Divergence-Free SPH (DFSPH)와 같은 방법들은 비압축성 유체에 존재하는 일정한 밀도 제약 조건과 발산 없는 속도 제약 조건을 모두 해결하기 위해 두 개의 솔버를 사용하며(Bender and Koschier/2015/Divergence-free smoothed particle hydrodynamics), 수렴을 개선할 수 있지만, 시스템에 기존 오차가 있을 때 높은 에너지 주입의 위험도 있습니다. 이러한 고에너지 상태로부터 복구하는 것은 충분한 반복 횟수나 충분히 낮은 타임스텝 없이는 불가능한 것으로 입증되었습니다.

PBF(Macklin and Müller/2013/Position based fluids) 및 Semi-Implicit SPH (SISPH)(He et al./2025/A Semi-Implicit SPH Method)와 같은 위치 기반 솔버는 큰 타임스텝과 적은 반복 횟수로 시뮬레이션을 안정적으로 유지하기 위한 유망한 접근 방식이며, 일반적으로 밀도 변동을 최소화할 수 있는 최적의 위치를 찾을 때 과도한 수정을 방지하기 위해 서브스텝핑이나 신중한 스텝 크기 계산을 활용합니다. 우리 방법과 유사하게, SISPH는 암시적 오일러의 변분 형태에서 파생되었지만, 시스템의 에너지를 최소화하기 위해 준-암시적 연속 치환 방법을 사용하는 반면, 우리의 솔루션은 VBD(Chen et al./2024/Vertex Block Descent)와 유사하게 지역적 뉴턴 반복을 사용하는 2차 하강 공식을 사용합니다.

SPH의 다른 개선 사항들은 고점성 유체 운동 처리(Peer et al./2015/An implicit viscosity formulation for SPH fluids) 및 SPH 유체 시뮬레이션의 인장 불안정성 개선(He et al./2020/A Variational Staggered Particle Framework; Liu et al./2024/A Dual-Particle Approach)과 같이 비압축성 유체에 대한 SPH 접근법의 다른 미해결 문제들을 다룹니다. 실제로, 우리의 연구는 이러한 개선 사항들과 병행하여 작동하며, 유체 운동을 더욱 개선하기 위해 우리 솔버와 함께 사용될 수 있습니다.

저희는 Gast와 Schroeder (2015)의 연구를 따라 암시적 오일러의 변분 형태를 사용하여 SPH에 대한 암시적 적분을 공식화합니다. 이 접근법은 강한(stiff) 제약 조건을 사용하여 비압축성을 시뮬레이션할 수 있게 해줍니다. 결과로 나오는 최소화 문제는 위치 업데이트와 함께 완화된 자코비 반복법(relaxed Jacobi iterations)을 사용하여 해결합니다. 마지막으로, 안정성을 유지하기 위해 2차 항에 대한 근사를 논의하고, 높은 강성도의 해에서 발생하는 에너지 주입을 상쇄하기 위한 감쇠 공식을 소개합니다.

이번 섹션에서는 먼저 변분 형태를 이용한 암시적 오일러 적분을 유도하고(섹션 3.1), 그 다음으로 임의로 큰 강성 값을 처리하기 위해 변분 에너지를 어떻게 조정하는지 설명합니다(섹션 3.2). 이어서 수치 계산 결과를 통합하고 결과 위치를 업데이트하는 방법을 보여줍니다(섹션 3.4). 마지막으로, 불확정 행렬로 인한 불안정성을 더욱 방지하기 위해 헤시안 근사 기법(섹션 3.5)과 높은 강성도의 해에서 주입되는 에너지를 상쇄하기 위한 인공 감쇠 공식(섹션 3.6)을 설명합니다.

암시적 오일러의 변분 형태는 각 타임 스텝에 대해 변분 에너지 $Ψ$를 최소화하는 것을 목표로 합니다.

암시적 오일러의 변분 형태는 각 타임 스텝마다 변분 에너지 Ψ를 최소화하는 것을 목표로 합니다. 위치 기반 공식에서는 이를 다음과 같이 쓸 수 있습니다.

모든 입자의 위치 $x$는 변분 에너지를 최소화함으로써 결정됩니다.

여기서 첫 번째 항은 운동량 포텐셜로, 타임 스텝 크기 $h$와 모든 입자의 위치 $x$와 그들의 관성 위치 $y$의 차이에 대한 질량 가중 노름 $||·||_M$을 포함합니다.

$y$는 타임 스텝 시작 시점의 입자 위치 $x_t$와 속도 $v_t$, 그리고 상수 가속도 항 $a^★$를 사용하여 계산됩니다.

마지막 항 $E(x)$는 총 내부 포텐셜이며, $a^★$에 포함되지 않은 모든 내력의 포텐셜 에너지를 포함합니다. 이는 각 입자 $i$ 주변에서 계산된 모든 지역적 내부 포텐셜 항 $E_i(x)$의 합으로 쓸 수 있습니다.

중력으로 인한 가속도 및 타임 스텝 동안 일정한 다른 외력들은 $a^★$에 포함되어야 합니다. 점성력이 타임 스텝 내에서 일정하게 유지된다고 가정함으로써, 우리는 이를 $a^★$에 포함시켜 $E(x)$에 남는 것은 압력으로 인한 포텐셜 에너지가 되도록 합니다.

일반적인 SPH 공식에 따라, 우리는 압력 에너지를 이차 포텐셜로 모델링할 수 있습니다.

여기서 $k$는 강성(stiffness)이고 $C_i(x)$는 밀도 제약으로, 이는 SPH에서 종종 상태 방정식이라고 불립니다. 이전 연구에는 상태 방정식에 대한 다양한 공식이 있지만, 우리는 다음과 같은 간단한 형태를 선택합니다.

여기서 $\rho_i$는 입자 $i$ 주변의 유체 밀도이고 $\rho$는 기본 밀도입니다. 이 간단한 제약 공식은 정지 밀도로부터의 편차에 페널티를 부과하며, 우리의 암시적 해법과 함께 실제적으로 잘 작동합니다. 우리는 이 공식이 (매우) 비선형적인 항을 사용하는 다른 형태들보다 더 나은 안정성을 제공한다는 것을 발견했습니다. 하지만, 이를 맹목적으로 적용해서는 안 됩니다. 특히, 유체 표면에서는 밀도가 자연스럽게 0으로 사라져야 합니다. 음압으로 인한 뭉침 현상을 피하기 위해, 일반적인 해결책은 제약을 0으로 클램핑하는 것이며(Koschier/2022/SPH Survey), 우리는 모든 테스트에서 이 접근법을 사용합니다.

수식 1의 최소화 문제를 풀기 위해, 우리는 다른 모든 입자가 고정되어 있다고 가정함으로써 버텍스 블록 하강(VBD) [Chen et al./2024/VBD]과 유사한 접근 방식을 따릅니다.

그러면 각 입자에 대해 전역 변분 에너지는 입자당 지역 변분 에너지로 축소될 수 있습니다.

여기서 $m_i$, $\mathbf{x}_i$, $\mathbf{y}_i$는 각각 입자 $i$의 질량, 위치, 관성 위치이며, $\mathcal{P}_i$는 SPH 커널 반경 내에 있는 입자 $i$ 주변의 모든 입자 집합입니다.

따라서 $\mathcal{P}_i$는 자신의 밀도를 평가하기 위해 입자 $i$를 고려하는 모든 입자를 포함합니다.

VBD에서와 마찬가지로, $\mathbf{x}_i$만을 움직여 지역 변분 에너지 $\Psi_i(\mathbf{x})$를 최소화하면 전역 변분 에너지 $\Psi(\mathbf{x})$의 감소가 보장됩니다.

이는 가우스-자이델 반복법을 통해 효율적인 수치 적분 방법을 형성합니다.

하지만 가우스-자이델은 일반적인 SPH 조건 하에서 $\mathcal{P}_i$가 수십만 개의 입자를 포함할 수 있어 병렬 계산의 잠재력을 심각하게 제한하므로 엄청나게 비용이 많이 들 수 있습니다.

더욱이, 각 입자를 이동시킨 후에는 모든 이웃 입자의 밀도를 업데이트해야 하는데, 이는 상당한 계산 비용을 발생시키고 동일한 이웃 입자가 업데이트될 때 잠재적인 데이터 해저드를 유발합니다.

이는 병렬화의 추가적인 감소나 성능을 저해하는 원자적 연산의 사용을 의미할 것입니다.

따라서 우리는 SPH 기반 방법에서 일반적인 완화된 자코비 기법을 선택합니다 [Ihmsen et al./2014/IISPH; Solenthaler and Pajarola/2009/Predictive-corrective SPH].

모든 입자에 대한 위치 업데이트 $\Delta \mathbf{x}_i$를 계산한 후, 병렬적으로 $\mathbf{x}_i \leftarrow \mathbf{x}_i + \Delta \mathbf{x}_i / 2$를 사용하여 그 절반을 적용합니다.

VBD와 유사하게, 우리는 입자 위치를 뉴턴 방법을 사용하여 조정하며, 이는 2차 선형 시스템을 푸는 것을 포함합니다.

여기서 $f_i$는 변분 형태의 음의 공간 도함수로부터 계산된 힘입니다.

그리고 $H_i$는 변분 형태의 헤시안입니다.

SPH는 스무딩 커널 $W$를 사용하여 이웃 입자들의 위치와 질량으로부터 입자 밀도 $\rho_i$를 계산합니다.

따라서, 밀도 제약의 도함수는 다음과 같이 커널 함수의 도함수를 사용하여 계산할 수 있습니다.

여기서 $W_{ij} = W(\mathbf{x}_i - \mathbf{x}_j)$입니다. 저희 구현에서는 Koschier et al./2019/SPH Techniques의 3차 스플라인 커널을 사용하지만, 저희 공식은 어떤 SPH 커널 공식과도 작동합니다.

위치 변화량은 3x3 행렬의 역행렬로 계산됩니다.

여기서 첫 번째 항은 항상 대칭 양의 정부호(symmetric positive-definite)이지만, 두 번째 항은 특히 입자 부족 문제(particle deficiency problem)를 해결하기 위해 일반적인 압력 클램핑(pressure clamping) 기법을 사용할 때 양의 정부호와 가역성을 보장할 수 없습니다.

저희는 헤시안을 근사하여 이 문제를 피하는데, 두 번째 항을 열(column)의 노름(norm)에 기반한 대각 근사(diagonal approximation)로 대체합니다 Andrews et al. 2017.

이러한 처리는 근사된 헤시안이 항상 가역적임을 보장합니다.

암시적 통합과 위치 기반 업데이트는 수렴을 방해하는 큰 타임스텝과 제한된 반복 횟수에도 불구하고 안정적인 SPH 시뮬레이션을 가능하게 합니다. 이전 타임스텝에서 남은 제약 조건 오차는 수치적 폭발로 이어지지는 않지만, 밀도 제약 조건의 강성(stiffness)에 따라 에너지를 주입할 수 있습니다. 암시적 오일러 방법은 수치적 감쇠를 도입하는 것으로 알려져 있지만, 매우 강한 제약 조건과 제한된 반복 횟수로 비압축성 유체를 시뮬레이션할 때 주입되는 에너지가 이 효과를 압도하여 시뮬레이션이 움직임 없이 안정된 해에 도달하는 것을 막을 수 있습니다. 이는 일부 응용 프로그램에서는 바람직하지 않을 수 있습니다. 따라서 저희는 더 낮은 강성을 사용하는 해와 비교하여 유체로부터 운동 에너지를 추출하는 인공 감쇠 공식을 포함합니다.

저희의 감쇠 공식은 운동 에너지를 제거함으로써 입자의 최종 속도만을 수정하며, 최소한의 영향을 미치도록 설계되었습니다. 감쇠는 애니메이션이 거의 멈추려고 할 때만 적용되며, 수치적 오류로 인해 도입된 여분의 운동 에너지를 제거하려고 시도합니다.

먼저 더 큰 컴플라이언스 파라미터 $\alpha^$를 사용하여 각 입자에 대한 대안 위치 $x_i^$를 계산하는 것으로 시작합니다 (테스트에서는 $\alpha^=1/1000$을 사용했습니다). 실제로는 솔버의 마지막 반복 동안, 마지막 반복 시작 시의 위치 $x$에서 시작하여 $x_i^$를 계산하기 위해 단일 반복을 사용합니다. 이런 방식으로 $x_i^$의 계산은 최소한의 오버헤드를 수반합니다.

마지막 반복이 끝날 때, 두 위치 $x_i$와 $x_i^$가 커널 반경 $r$의 사용자 지정 배수 $\beta$보다 가까울 경우에만 감쇠를 적용합니다 (테스트에서는 $\beta=60$을 사용했습니다). 이는 높은 강성을 통해 많은 양의 운동량이 주입되었음을 의미합니다. 이 경우, 각 입자에 대해 최종 위치와 대안 위치를 사용하여 두 개의 속도를 다음과 같이 계산합니다.

그런 다음 입자의 운동 에너지 $K_i$를 이 대안 속도를 사용한 대안 운동 에너지 $K_i^$와 비교합니다. 만약 $K_i \le K_i^$이면, 어떠한 감쇠도 적용하지 않습니다. 그렇지 않으면, 입자의 운동 에너지를 대안에 더 가깝게 만듭니다. 따라서 최종 속도는 다음과 같이 계산됩니다.

여기서 감쇠 가중치 $d$는 다음을 사용하여 계산됩니다.

입력: $x_t$: 타임스텝 시작 시의 위치.

입력: $v_t$: 타임스텝 시작 시의 속도.

입력: $a^\star$: 타임스텝 시작 시의 비압력 힘.

출력: $x_{t+1}, v_{t+1}$, 타임스텝 종료 시의 위치와 속도

모든 입자에 대해

1: $y \leftarrow x_t + h v_t + h^2 a^\star$

2: $x_{\text{guess}} \leftarrow y$

3: $l = 0$

4: while $l < MAX_{\text{iterations}}$ do

5:     $x_{\text{guess}}$를 기반으로 모든 밀도 및 밀도 그래디언트 업데이트

6:     for 모든 입자 $i$ do

7:         수식 10을 사용하여 $f_i$ 계산

8:         수식 11을 사용하여 $H_i$ 계산

9:         $\Delta x_i \leftarrow H_i^{-1} f_i$

10:         $x_{\text{new}_i} \leftarrow x_{\text{guess}_i} + \Delta x_i / 2$

11:     end for

12:     $x_{\text{guess}} \leftarrow x_{\text{new}}$

13:     $l = l + 1$

14: end while

15: $x_{t+1} \leftarrow x_{\text{guess}}$

16: 수식 17을 사용하여 $v_{t+1_i}$ 설정

저희는 NVIDIA RTX 4090에서 CUDA를 사용하여 저희의 방법(알고리즘 1 참조)을 구현했습니다. 별도로 명시하지 않는 한, $\alpha=0$과 저희의 감쇠를 사용합니다. 입자 크기는 직경 0.5, 커널 크기는 1로 설정됩니다. 저희는 IPBF 방법을 PBF [Macklin and Müller 2013], IISPH [Ihmsen et al. 2014], SISPH [He et al. 2025], 그리고 DFSPH [Bender and Koschier 2015]와 비교합니다. 모든 방법은 반복당 계산 비용이 비슷하지만, 한 방법이 유체 밀도를 유지하지 못하고 입자들이 더 많은 이웃을 갖기 시작하면 계산 비용이 증가할 수 있습니다. 따라서 저희는 다른 방법을 비교할 때 단순히 동일한 반복 횟수를 사용하기로 했습니다. 한 가지 예외는 DFSPH로, 최소 1회의 발산 반복과 2회의 밀도 반복이 필요합니다(반복 횟수를 보고할 때 1+2로 표기).

그림 2와 그림 3은 두 유체 기둥이 중앙에서 충돌하여 얇은 제트를 형성하는 Double Dam Break 예제의 결과를 보여줍니다. 그림 2에서 보듯이, 모든 방법이 4회 반복(DFSPH는 2+2)으로 안정적일 때 비슷한 결과를 생성합니다. 이 경우, DFSPH는 약간의 수치적 불안정성을 보이지만 밀도 보존은 탁월하며, SISPH는 밀도 변동으로 인해 약간의 추가 에너지 주입이 있고, IISPH는 추가 압축으로 인해 에너지를 소산하며, PBF는 다른 방법들과 비교하여 제트의 높이에서 볼 수 있듯이 가장 많은 에너지를 소산합니다.

우리는 알고리즘 1에 보이는 바와 같이 NVIDIA RTX 4090에서 CUDA를 사용하여 우리의 방법을 구현했습니다. 별도로 명시하지 않는 한, 우리는 $𝛼= 0$과 우리의 댐핑을 사용합니다. 입자 크기는 직경 0.5, 커널 크기는 1로 설정했습니다. 우리는 우리의 IPBF 방법을 PBF Macklin/2013/Position based fluids, IISPH Ihmsen/2014/Implicit Incompressible SPH, SISPH He/2025/A Semi-Implicit SPH Method, 그리고 DFSPH Bender/2015/Divergence-free smoothed particle hydrodynamics와 비교합니다. 모든 방법은 반복당 계산 비용이 비슷하지만, 유체 밀도를 유지하지 못하고 입자들이 더 많은 이웃을 갖게 되면 계산 비용이 증가할 수 있습니다. 따라서, 우리는 다른 방법들을 비교할 때 단순히 동일한 반복 횟수를 사용하기로 했습니다. 한 가지 예외는 DFSPH로, 최소 1회의 발산 반복과 2회의 밀도 반복이 필요합니다 (반복 횟수를 보고할 때 1+2로 표기됨).

그림 2와 그림 3은 우리의 Double Dam Break 예제의 결과를 보여주는데, 여기서 두 개의 유체 기둥이 중앙에서 충돌하여 얇은 제트를 형성합니다. 그림 2에서 볼 수 있듯이, 모든 방법은 4회 반복(DFSPH의 경우 2+2)으로 안정적일 때 비교 가능한 결과를 생성합니다. 이 경우, DFSPH는 약간의 수치적 불안정성을 보이지만 뛰어난 밀도 보존을 나타내고, SISPH는 밀도 변동으로 인한 약간의 추가 에너지 주입이 있으며, IISPH는 과도한 압축으로 인해 에너지를 소산시키고, PBF는 가장 많은 에너지를 소산시키는데, 이는 다른 것들과 비교하여 제트의 높이에서 볼 수 있습니다.

그림 2의 매 프레임 시작 시점에서의 밀도 오차를 측정하면, 그림 4에 그려진 바와 같이 모든 방법에서 안정적인 애니메이션 결과를 더욱 명확하게 보여줍니다. 보충 비디오의 애니메이션에서 예상할 수 있듯이, 그래프는 DFSPH가 안정적일 때 우수한 밀도 보존을 보이며, 심지어 우리의 방법을 능가하는 것을 보여줍니다. 반면 SISPH와 IISPH는 더 높은 밀도 오차와 변동을 보입니다. PBF는 가장 큰 밀도 오차와 변동을 가지며, 이는 애니메이션이 정지 상태에 가까워질 때조차도 밀도 오차의 큰 변화로 나타납니다.

그림 3은 시간 스텝당 2회 반복(DFSPH의 경우 1+2)의 결과를 보여줍니다. 이 경우, DFSPH와 IISPH는 수치적 폭발로 인해 실패하지만, 모든 위치 기반 방법들은 이전과 유사한 결과를 생성합니다.

안정성. 그림 5는 안정성에 대한 스트레스 테스트를 보여줍니다. 이 예제에서, 우리는 초기 상태가 기준 밀도의 7배인 유체 구를 초기화하고 동일한 시간 스텝 크기와 동일한 반복 횟수로 시뮬레이션합니다. 단, DFSPH는 발산 해결을 위해 +3회의 반복이 추가로 주어집니다. 그럼에도 불구하고, 다른 방법들과 달리 DFSPH는 큰 밀도 편차에서 회복하지 못하고, 시뮬레이션 시작 직후 수치적 폭발을 경험합니다. IISPH는 안정성을 유지하지만, 기준 밀도로 회복할 수 없기 때문에 상당한 부피 손실을 대가로 합니다. 우리와 유사하게, PBF와 SISPH는 상당한 부피 손실 없이 회복할 수 있지만, 일관된 흐름을 구성하는 데 실패하여 많은 개별 유체 입자들이 날아다닙니다. 대조적으로, 우리 방법은 부피 손실이나 비일관적인 흐름 없이 이 상당한 밀도 편차로부터 회복합니다. 부피 손실은 종종 안정성과의 트레이드오프 관계에 있기 때문에 이는 주목할 만합니다. 더 강성인 시스템은 압축을 줄일 수 있지만, 일반적으로 수치적으로 덜 안정적입니다. 우리 방법은 둘 다 달성할 수 있습니다. 사실, 우리 방법은 모든 테스트에서 무조건적인 안정성을 보입니다.

우리는 그림 6에서 우리 방법의 이러한 동작을 보여주는데, 여기서 우리는 그림 2의 동일한 Double Dam Break 예제를 사용하여 시간 스텝 크기 $h$를 점진적으로 증가시키면서 시간 스텝당 2회의 반복만 사용합니다. 우리 방법은 $h$를 1/240에서 1/60초로 높여도 거의 변화를 보이지 않습니다. 심지어 1/30초의 시간 스텝 크기에서도, 우리 방법은 눈에 띄는 압축을 대가로 안정성을 유지합니다. 이 테스트는 우리 방법이 많은 프레임에 걸쳐 축적된 상당한 잔여 오차를 가지고 작동하며, 큰 시간 스텝 크기에 대해 단 2회의 반복으로는 밀도 오차를 수정할 수 없지만, 여전히 안정성을 유지함을 보여줍니다. 일반적으로, 우리 방법은 필요할 때 비압축성을 안정성과 맞바꾸며, 따라서 극한의 시간 스텝 크기에서도 안정성을 유지합니다.

비압축성. 이러한 트레이드오프가 필요하지 않을 때, 우리 방법은 자연스럽게 우수한 비압축성을 나타냅니다. 우리는 이를 2D로 시뮬레이션된 큰 정적 물기둥을 사용하여 그림 7에서 보여줍니다. 동일한 계산 시간 하에서, 우리 방법은 DFSPH와 유사한 압축률을 보이지만, PBF, SISPH, IISPH는 상당히 압축됩니다. 그림 8에서는 모든 방법에 대해 처음 500 프레임 동안의 평균 밀도를 보여줍니다. 우리 IPBF 방법은 빠르게 정상 상태에 도달하는 반면, 다른 방법들은 더 높은 압축 외에도 더 많은 진동을 보인다는 점에 주목하십시오.

동일한 동작이 3D에서도 관찰될 수 있습니다. 그림 11에서는 물웅덩이 위로 물 블록을 떨어뜨립니다. 이는 접촉 시 많은 양의 압축 오차를 발생시킵니다. 그림 12는 이 접촉 직후의 시뮬레이션 상태를 다른 솔버들과 비교하여 보여줍니다. 여기서, 비압축성을 강제하려는 시도에서 DFSPH는 안정성을 유지하지 못하고 접촉 시 폭발합니다. PBF와 IISPH는 안정성을 유지했지만, 과도한 압축을 보였고, 이는 또한 높은 압축 상태에서는 각 입자가 주변의 커널 반경 내에서 훨씬 더 많은 입자들과 상호작용하기 때문에 계산 비용의 상당한 증가로 이어집니다. 이 예제에서 SISPH도 안정성을 유지하지만, 큰 오차 잔여물로 인해 더 큰 에너지 주입을 보입니다.

헤시안 근사는 수치적 안정성을 달성하는 데 매우 중요합니다. 그림 15에서 볼 수 있듯이, 우리의 헤시안 근사 없이는 시뮬레이션이 비교적 간단한 경우에도 종종 빠르게 불안정해집니다.

우리는 큰 타임스텝에서도 뛰어난 안정성을 보이면서 낮은 밀도 오차를 유지하는, SPH 유체 해결을 위한 암시적 위치 기반 유체 방법을 소개했습니다. 우리는 우리 방법이 이전 연구들을 능가할 수 있음을 보여주며, 비슷한 시뮬레이션 품질을 보이면서 극단적인 타임스텝에서도 안정적으로 유지됩니다.

안정성을 유지함에도 불구하고, 우리 방법은 비압축성 유체에 존재하는 높은 강성으로 인한 에너지 주입 때문에 아티팩트를 보일 수 있습니다. 이는 주로 $𝛼$값이 0일 때 발생하는 잔여 오차 때문이며, 이로 인해 우리 방법은 수렴 보장을 제공하지 못합니다. 향후 연구에서는 시뮬레이션을 위해 사용자가 직접 튜닝할 필요가 없는 대안적인 감쇠 공식을 제시하여 이 문제를 해결할 수 있습니다.

Fig 1: NVIDIA RTX 4090 GPU에서 저희 방법을 사용하여 구 위로 떨어지는 백만 개의 유체 입자를 시뮬레이션한 모습입니다. 시뮬레이션은 프레임당 평균 159밀리초로 실행되었습니다. 비압축성 유체의 효율적인 시뮬레이션은 여전히 어렵고 미해결된 문제입니다. 이전 연구들은 종종 비압축성, 안정성, 비용 사이에서 다양한 절충을 합니다. 그러나 세 가지 모두를 얻는 경우는 드뭅니다. 본 논문에서는 비압축성에 접근하기 위해 특별히 공식화된 변분 에너지를 최적화하기 위해 2차 암시적 하강 기법을 사용하는 새로운 비압축성 Smoothed Particle Hydrodynamics (SPH) 기법을 소개합니다. 저희는 저희 방법이 최소한의 계산 비용으로 비압축성과 안정성 모두에서 우수함을 보여줍니다. 또한, 저희 방법은 극한의 타임 스텝에서도 무조건적으로 안정적이어서 대화형 애플리케이션에 적합함을 보여줍니다. smoothed particle hydrodynamics (SPH)에 기반한 라그랑지안 입자는 컴퓨터 그래픽스에서 유체 시뮬레이션을 위한 인기 있는 기술이었습니다. 비압축성 유체를 시뮬레이션하는 데 자주 사용되지만, SPH로 비압축성을 달성하는 것은 근본적인 과제입니다. 이는 SPH 공식에서 비압축성이 유체 밀도를 유지하는 강한 제약 조건에 해당하기 때문입니다. 결과적으로, 명시적 속도 적분은 안정성을 위해 매우 작은 타임 스텝을 필요로 하며, 고도로 비선형적인 제약 조건을 사용하여 비압축성을 달성하기 위해서는 훨씬 더 작은 타임 스텝이 필요합니다. 따라서 대부분의 최근 연구는 입자별로 제약 강성을 적응적으로 조정하는 준-암시적 속도 적분의 변형을 채택했습니다. 불행히도 이러한 접근 방식은 복잡한 시뮬레이션에서 안정성 문제나 과도한 압축으로 여전히 어려움을 겪으며 비교적 작은 타임 스텝을 요구합니다. 이는 병렬 컴퓨팅을 통한 계산 성능이라는 SPH 기반 방법의 가장 큰 장점 중 하나를 저해합니다. 한 가지 주목할 만한 예외는 position-based fluids (PBF) 방법으로, SPH를 위치 기반 동역학과 유사한 제약 공식으로 변환하여 더 큰 타임 스텝에서 뛰어난 안정성을 제공합니다.

Fig 2: 타임 스텝당 4번의 이터레이션을 사용한 Double Dam Break 예제. PBF, SISPH, IISPH, DFSPH(2+2 이터레이션), 그리고 저희 방법으로 실행되었습니다. 이 그림의 모든 방법은 h=1/240초를 사용하여 달성된 60ms의 계산 예산으로 실행됩니다. 표시된 방법은 PBF, IISPH, SISPH, DFSPH, IPBF (저희 방법)입니다.

Fig 2: 이 경우, DFSPH는 약간의 수치적 불안정성을 보이지만 밀도 보존은 우수하며, SISPH는 밀도 변동으로 인해 약간의 추가 에너지가 주입되고, IISPH는 추가 압축으로 인해 에너지를 소산하며, PBF는 다른 방법들과 비교하여 제트의 높이에서 볼 수 있듯이 가장 많은 에너지를 소산합니다. 모든 방법에서 안정적인 애니메이션 결과로부터 추가적인 하이라이트를 얻기 위해 매 프레임 시작 시 그림 2의 밀도 오차를 측정했습니다. 그래프는 PBF, IISPH, SISPH, IPBF (저희 방법), DFSPH에 대한 평균 상대 오차 대 프레임을 보여줍니다.

Fig 3: 타임 스텝당 2번의 이터레이션을 사용한 Double Dam Break 예제. 단, DFSPH는 추가적인 발산 이터레이션을 사용하여(1+2 이터레이션) 실행되었습니다. 이 경우, IISPH와 DFSPH는 수치적 안정성을 유지하지 못합니다. 저희는 CUDA를 사용하여 NVIDIA RTX 4090에서 저희 방법을 구현했습니다. 별도로 명시되지 않는 한, α=0과 저희의 댐핑을 사용합니다. 입자 크기는 직경 0.5, 커널 크기는 1로 설정되었습니다. 저희는 IPBF 방법을 PBF, IISPH, SISPH, DFSPH와 비교합니다. 모든 방법은 이터레이션당 계산 비용이 비슷하지만, 방법이 유체 밀도를 유지하지 못하고 입자들이 더 많은 이웃을 갖기 시작하면 계산 비용이 증가할 수 있습니다. 따라서 저희는 다른 방법을 비교할 때 단순히 동일한 이터레이션 횟수를 사용하기로 했습니다. 한 가지 예외는 DFSPH로, 최소 1번의 발산 이터레이션과 2번의 밀도 이터레이션이 필요하며, 이터레이션 횟수를 보고할 때 1+2로 표기합니다.

Fig 4: 그림 2의 Double Dam Break 예제에 대한 여러 방법의 밀도 평균 상대 오차를 그림 4에 나타냈습니다. 보충 비디오의 애니메이션에서 예상할 수 있듯이, 그래프는 DFSPH가 안정적일 때 저희 방법을 능가할 정도로 우수한 밀도 보존을 보이며, SISPH와 IISPH는 더 높은 밀도 오차와 변동을 보입니다. PBF는 가장 큰 밀도 오차와 변동을 가지며, 이는 애니메이션이 정지 상태에 가까워질 때도 밀도 오차의 큰 변화로 나타납니다.

Fig 2: 저희 방법은 h를 1/240초에서 1/60초로 높여도 거의 변화를 보이지 않습니다. 타임 스텝 크기가 1/30초일 때조차도, 저희 방법은 눈에 띄는 압축을 대가로 안정성을 유지합니다. 이 테스트는 저희 방법이 많은 프레임에 걸쳐 축적된 상당한 잔차로 작동하며, 큰 타임 스텝에 대해 단 2번의 이터레이션으로는 밀도 오차를 수정할 수 없지만 여전히 안정적으로 유지됨을 보여줍니다. 일반적으로, 저희 방법은 필요할 때 안정성을 위해 비압축성을 희생하므로 극한의 타임 스텝 크기에서도 안정성을 유지합니다. 비압축성에 관해서, 이러한 절충이 필요하지 않을 때, 저희 방법은 자연스럽게 우수한 비압축성을 나타냅니다. 이는 2D로 시뮬레이션된 큰 정적 물기둥을 사용한 그림 7에서 보여줍니다.

Fig 5: 정지 밀도의 7배로 중앙에 압축된 유체 공으로 시작하는 압축 안정성 테스트입니다. 저희 IPBF 방법은 안정된 상태로 회복하는 반면, DFSPH는 폭발하고, IISPH는 높은 밀도 오차에서 회복하지 못하며, PBF와 SISPH 모두 불안정한 입자들을 남깁니다. 모든 예제는 h=1/120초와 3번의 이터레이션(DFSPH는 3+3)으로 실행되었습니다. 상단 행은 애니메이션의 초기 프레임을 보여줍니다. 이미지들은 h=1/240초, 1/120초, 1/60초, 1/30초에 대한 결과를 보여줍니다.

Fig 6: 저희 방법은 타임 스텝 크기를 h=1/30초까지 늘리고 이터레이션 횟수를 타임 스텝당 2회로 제한하더라도 더블 댐 붕괴에 대해 안정적인 결과를 생성합니다. 큰 타임 스텝 크기와 제한된 이터레이션으로 인한 오차는 과도한 압축으로 이어집니다.

Fig 7: 2D에서 물기둥을 떨어뜨려, 각 방법이 8번의 이터레이션(DFSPH는 4+4)과 α=10⁻⁶(SISPH는 μ=10⁶)에서 얼마나 많은 압축을 보이는지 측정했습니다. 수렴된 상태에서, IPBF는 시각적으로 가장 적은 압축을 보입니다. 그래프는 PBF, IISPH, SISPH, DFSPH, IPBF (저희 방법)에 대한 평균 상대 밀도 대 프레임을 보여줍니다.

Fig 8: 이 그래프는 그림 7과 보충 비디오에 나온 2D 압축 테스트에서 애니메이션 과정 동안의 입자 평균 밀도를 보여줍니다. IPBF는 매우 빠른 수렴과 정지 밀도(1)로부터의 최소한의 편차를 보입니다. 동일한 계산 시간 하에서, 저희 방법은 DFSPH와 유사한 압축률을 보이며, 반면 PBF, SISPH, IISPH는 상당히 압축됩니다. 그림 8에서는 모든 방법에 대해 처음 500 프레임 동안의 평균 밀도를 보여줍니다. 저희 IPBF 방법은 빠르게 안정 상태에 도달하는 반면, 다른 방법들은 더 높은 압축 외에도 더 많은 진동을 보입니다. 동일한 현상이 3D에서도 관찰됩니다. 그림 11에서는 물웅덩이 위로 물 블록을 떨어뜨립니다. 이는 접촉 시 많은 양의 압축 오차를 발생시킵니다. 그림 12는 이 접촉 직후 다른 솔버들을 비교한 시뮬레이션 상태를 보여줍니다. 여기서 비압축성을 강제하려는 시도에서 DFSPH는 안정성을 유지하지 못하고 접촉 시 폭발합니다. PBF와 IISPH는 안정성을 유지했지만 과도한 압축을 보였고, 이는 또한 계산 비용의 상당한 증가로 이어집니다. 왜냐하면 높은 압축 상태에서는 각 입자가 커널 반경 내에서 훨씬 더 많은 입자들과 상호작용하기 때문입니다. 이 예제에서 SISPH도 안정적으로 유지되지만, 큰 오차 잔차로 인해 더 큰 에너지 주입을 보입니다.

Fig 9: 3D 수렴 테스트를 실행하는 데 사용된 프레임입니다. 큰 물 블록이 아래의 물웅덩이에 떨어지는 순간입니다. 그래프는 PBF, IISPH, μ=10⁸인 SISPH, DFSPH, α=10⁻⁸인 IPBF, α=0인 IPBF에 대한 평균 상대 오차 대 이터레이션을 보여줍니다.

Fig 10: 저희 방법은 그림 9에 표시된 동일한 초기 상태에서 시작할 때 더 적은 이터레이션으로 더 낮은 밀도 오차로 수렴합니다. 그림 10에서 수렴 테스트를 제시하며, 증가하는 이터레이션 횟수에 따른 한 타임 스텝 동안의 유체 밀도 평균 상대 오차를 보여줍니다. 이 테스트의 모든 방법은 그림 9에 표시된 정확한 시뮬레이션 상태에서 시작하며, 이는 그림 11의 애니메이션에서 물 블록이 웅덩이와 충돌하는 프레임에 해당합니다. 그림 10에서 저희 IPBF 방법이 이 테스트에서 10회 미만의 이터레이션으로 낮은 평균 상대 오차로 빠르게 수렴하는 것을 주목하십시오. 0이 아닌 알파 파라미터를 사용하면 변분 에너지의 운동량 항을 고려하기 위해 약간 더 큰 밀도 오차가 관찰됩니다. DFSPH는 속도 발산 해결을 위해 2번의 추가 이터레이션을 포함하며, 이는 수렴을 지연시킵니다. 그렇지 않으면 이 방법도 빠르게 수렴하지만 더 큰 오차로 수렴합니다. 다른 방법들은 훨씬 느린 수렴을 보입니다. 이 그래프에는 표시되지 않았지만, SISPH는 추가 이터레이션을 통해 DFSPH의 오차 수준 이하로 오차를 계속 줄이지만, 1000번의 이터레이션 후에도 저희 IPBF 방법보다 높은 오차를 유지합니다. 속도에 관해서, 저희 방법은 이터레이션당 계산 비용이 비슷하면서도 안정적인 시뮬레이션을 생성하기 위해 더 큰 타임 스텝 크기를 허용하므로, 저희 방법은 이전 연구보다 훨씬 빠를 수 있습니다.

Table 1: 유사한 평균 밀도 오차에서의 성능 비교. 각기 다른 방법들에 대해 스텝당 다른 이터레이션 횟수와 프레임당 다른 스텝 크기를 사용하여 달성되었습니다. 프레임당 총 계산 시간은 NVIDIA RTX 4090 GPU에서 측정되었습니다.

Fig 11: 저희 IPBF 방법을 사용하여 시뮬레이션된 25만 개 입자의 Block Flop 예제입니다. 유체 블록이 웅덩이 위로 떨어집니다. 결과 시뮬레이션은 높은 밀도 오차 기간을 초래하지만, 2번의 이터레이션과 h=1/240의 타임 스텝 크기로 우리는 비압축성을 유지하며 안정성을 유지하고, 프레임당 평균 60ms로 실행됩니다. 표시된 방법은 PBF, IISPH, SISPH, DFSPH, IPBF (저희 방법)입니다.

Fig 12: 그림 11의 Block Flop 애니메이션에서 블록이 웅덩이에 부딪힌 직후의 프레임으로, 다른 방법들을 사용하여 시뮬레이션되었습니다. 모든 방법은 h=1/240초와 2번의 이터레이션을 사용하여 프레임당 60ms의 유사한 계산 시간을 가집니다. 단, DFSPH는 추가 계산 비용을 유발하는 2번의 추가 발산 이터레이션을 가졌습니다. 그럼에도 불구하고, DFSPH는 이 예제에서 완전히 불안정합니다. SISPH는 밀도 오차로 인해 에너지를 주입합니다.

Fig 13: α=0을 사용한 해결책과 저희의 인공 감쇠 공식을 사용한 동일한 해결책으로 작은 댐 붕괴를 비교합니다. 감쇠가 없으면 유체 표면이 더 거칠고 정지하지 않지만, 저희의 감쇠를 사용하면 일정 시간 후에 정지 상태로 안정됩니다. 보충 비디오에서 볼 수 있듯이, 저희의 감쇠 공식은 전체적인 유체 운동에 최소한의 영향을 미치며, 입자들이 거의 정지 상태에 가까워질 때만 운동을 감쇠시킵니다.

Fig 14: 다른 유체 압축성을 시뮬레이션하기 위해 다른 컴플라이언스 파라미터 α를 사용하여 저희 IPBF 방법으로 댐 붕괴를 시뮬레이션한 모습입니다. 이미지들은 헤시안 근사 사용 여부에 따른 결과를 비교합니다.

Fig 15: 저희의 헤시안 근사 없이는 수치적 안정성을 유지하지 못하는 댐 붕괴 시뮬레이션입니다. 컴플라이언스에 관하여, 저희 IPBF 방법은 다른 압축성 거동을 가진 유체를 시뮬레이션하기 위해 다른 컴플라이언스 파라미터 α와 함께 사용될 수 있습니다. 이는 그림 14에서 보여줍니다. 헤시안 근사에 관하여, 저희의 헤시안 근사는 수치적 안정성을 달성하는 데 중요합니다. 그림 15에서 볼 수 있듯이, 저희의 헤시안 근사 없이는 비교적 간단한 경우에도 시뮬레이션이 종종 빠르게 불안정해집니다.