물 표면 시뮬레이션을 위한 보존적 USCIP법

전, 세종; 송, 오영

doi:10.15701/kcgs.2019.25.5.21

J Korea Comput Graph Soc 2019; 25(5):21-30

pISSN: 1975-7883, eISSN: 2383-529X

DOI: https://doi.org/10.15701/kcgs.2019.25.5.21

Article

물 표면 시뮬레이션을 위한 보존적 USCIP법

전세종¹

, 송오영²^,^*

A Conservative USCIP Simulation Method for Shallow Water

Sejong Jeon¹

, Oh-young Song²^,^*

Author Information & Copyright ▼

¹세종대학교 소프트웨어학과

²세종대학교 소프트웨어학과

¹Department of Software, Sejong University, ypland@gmail.com

²Department of Software, Sejong University, oysong@sejong.ac.kr

^*corresponding author: Oh-young Song/Sejong University(oysong@sejong.ac.kr)

© Copyright 2019 Korea Computer Graphics Society. This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

Received: Nov 15, 2019; Revised: Nov 25, 2019; Accepted: Nov 29, 2019

Published Online: Dec 31, 2019

요약

이 논문은 물 표면을 효율적, 효과적으로 표현하기 위한 물리적 시뮬레이션 방법을 제안한다. 이 논문에서 표현하고자 하는 물은 깊이에 비해 너비가 매우 크고 상하 유동이 적은 상태로서, 이를 효율적으로 계산하기 위해 Navier-Stokes 방정식을 간략화한 천수방정식(shallow water equation)을 지배방정식으로 사용한다. 천수방정식의 대류항을 수치적으로 계산하기 위한 방법으로 기존의 Constrained Interpolation Profile(CIP)법을 개선하여, 수치적인 정확성을 높이고 물리량을 보존할 수 있는 Conservative Unsplit Semi-lagrangian CIP(CUSCIP)을 소개한다. 이 방법은 Kim 등이 제안한 USCIP[9]기법에서 사용하는 제약 조건에 적분값을 반영한 항을 추가하여 대류항을 계산한다. 실험결과를 통해, CUSCIP방법은 수치 소산(numerical dissipation)으로 인한 물리량 손실에 강건하여, 물결의 세밀함과 더불어 물결의 지속성이 향상됨을 알 수 있다.

Abstract

We propose a physical simulation method based on the shallow water equation(SWE) to represent water surface effectively. In this paper, the water which can be represented has a much larger width compared to the depth does not have a large vertical direction flow. In order to calculate the water flow efficiently, we start with the shallow water equation as the governing equation, which is a simplified version of the Navier-Stokes equation. In order to numerically calculate the advection term of the SWE, we introduce a new conservtive USCIP(CUSCIP) method which improves the Constrained Interpolation Profile (CIP) method to preserve the physical quantity while increasing the numerical accuracy. The proposed method is based on Kim et. al.’s Unsplit Semi-lagrangian CIP[9], and calculates advection term with additional constraints on term that consider integral values. The experimental results show that the CUSCIP method is robust to the loss of physical quantity due to numerical dissipation, which improves wave detail and persistence.

Keywords: 유체 시뮬레이션; 천수방정식; CIP법; 물리기반 애니메이션

Keywords: fluid simulation; shallow water equation; CIP; physical based animation

1. 서른

지구상에서 관찰할 수 있는 자연 현상들의 사실적 표현은 컴퓨터 그래픽스 분야의 필수적인 목표이다. 그중에서도 물은 자연에서 가장 흔하게 관찰할 수 있는 물질 중 하나로서, 이를 사실적으로 표현하기 위한 기술은 매우 중요한 과제로 연구되어왔다. 현재는 전산유체역학의 이론을 이용한 시뮬레이션 기법으로 이와 같은 물의 움직임을 보다 사실적으로 재현하고 있다.

유체를 시뮬레이션하기 위해서는 유체의 움직임을 의미하는 지배방정식을 풀어야 한다. 이는 Navier-Stokes 방정식이라 불리는 비선형 편미분 방정식으로 표현된다. 유체 시뮬레이션은 실제의 물리량을 이산화하고 그 값을 컴퓨터의 반복 연산으로 풀어내는 방식으로 유체의 움직임을 재현한다. 물이 깊이에 비해 너비가 훨씬 크고 상하 유동이 적은 상태일 경우 3차원의 Navier-Stokes 방정식은 2차원으로 간략화할 수 있으며 이를 천수방정식(Shallow Water Equation; SWE)이라 한다.

이 논문에서는 천수방정식의 대류항을 푸는 과정에서 발생하는 수치적 손실을 줄이기 위해 Constrained Interpolation Profile(CIP)법을 도입할 것이며, 천수에서 특히 중요한 물리량인 높이값을 보존하기 위해 일반적인 CIP를 개선한 방법, CUSCIP(Conservative Unsplit Semi-lagrangian CIP)를 제안한다.

2. 관련 연구

Kass와 Miller[3]는 2차원 천수방정식으로 높이 필드를 구하여 수면을 표현하였다. 이후 Layton과 Panne[13]은 안정성을 보장하는 시뮬레이션 방법을 개발했다. Yuksel, House, Keyser[14]는 접근 방식을 바꾸어 높이값을 변화시키는 파동 자체를 하나의 입자(particle)로 보고 그 입자를 시뮬레이션하는 방법을 사용하였다. Jeschke와 Wojtan[15]은 입자를 파동묶음(wave paket)으로 확장하여 다양한 주파수의 수면파를 시뮬레이션 하는 방법을 개발하였다. Lee와 Han[16]은 2차원 상에 체적값을 갖는 입자들을 SPH[17]로 시뮬레이션하여 천수방정식을 푸는 방법을 제안하였다. 수면파를 이산 푸리에 변환으로 해석하여 움직임을 표현하는 방법들도 제안되었다[18,19]. 한편, 쇄파나 분무 현상을 표현 못하는 천수방정식의 한계를 극복하기 위해 다양한 논문들이 발표되었다. O'Brien과 Hodgins[20], Thürey와 동료들[21]은 이 문제를 입자 시스템(particle system)을 도입하여 해결하였다. 즉, 높이 필드에서 쇄파나 분무 현상이 일어날 곳을 찾아 그 위치에 입자를 생성하는 방식이다. Irving과 동료들[22]은 쇄파나 분무 같은 세밀한 현상이 물 표면에서만 일어난다는 사실에 착안하여 2차원 격자 기반의 높이 필드와 일반적인 3차원 유체 시뮬레이션을 결합하는 방법을 제안하였다.

한편, 유체 시뮬레이션 대류항을 계산하는 가장 대표적인 수치 해법은 세미-라그랑지안 기법(Semi-Lagrangian method)이다. 전산 유체역학 분야에서 Courant, Isaacson, Rees[1]에 의해 제안되었고 그래픽스 분야에서는 Stam[2]이 최초로 소개하였다. 대류항을 수치적으로 풀 때 필요한 보간 과정에서의 정확도를 높이기 위해 Yabe와 동료들[4,5,6,7]은 CIP법을 고안하였다. CIP법은 물리량을 보간할 때 미분값을 함께 고려하는 방법이다. 최초의 CIP는 안정성을 보장하지 못하는 문제를 가지고 있었다. Song과 Shin, Ko는 CIP의 안정성 문제를 해결하기 위해 미분값을 보정하는 Monotonic CIP(MCIP)법을 제안하였다[8]. 이 방법은 Courant-Friedrichs-Lewy(CFL)조건에 무관하게 안정적이다. 이후 Kim 등은 [8]에서 야기되는 차원 분할 과정을 생략하기 위해 Unsplit Semi-lagrangian CIP(USCIP)법을 제안하였다[9]. 이 방법은 MCIP보다 정확도가 향상 되었으며, 계산 비용을 대폭 절감하였다. 물리량을 보존하기 위해 개발된 CIP-Conservative Semi-Lagrangian(CIP-CSL)법[10,11]은 적분값을 함께 고려한다. 2009년 Yabe와 Ogata는 이를 천수방정식에 적용하는 방법을 소개하였다[12].

3. 천수방정식의 수치해법

3.1 천수방정식

천수방정식을 유도하기 위해서는 다음과 같은 조건들이 필요하다.

유체는 얕은(shallow) 상태이다. 이는 수평 방향의 크기가 수직 방향의 크기보다 월등히 크다는 것을 의미한다.
유체는 유동이 적은, 안정적인 상태이다. 유체는 특정한 방향으로 흐르고 있지 않으며, 기본적으로 작용하는 외력은 중력뿐이다. 따라서 압력은 더 이상 미지의 값이 아니며, 깊이에 따라 선형으로 비례한다.
유체는 높이값을 갖는 자유표면으로 표현된다(height field). 따라서 쇄파(breaking wave)나 분무(spray)현상은 표현할 수 없다.

이와 같은 조건으로부터 유체는 수직 방향으로의 속도 변화를 무시할 수 있으며, 연직 방향에 대하여 정역학적(hydrostatic)상태라는 가정이 가능하다. 따라서 유체의 압력은 중력의 영향만을 고려하여 유체의 깊이에 선형으로 비례한다. 여기서 주의할 점은 수직 방향의 속도가 반드시 0은 아니라는 점이다. 속도의 수직 성분은 유체가 존재하는 지면의 높이와 중력에 의해 발생하며, 수평 방향의 속도가 깊이 변화에 영향을 받지 않을 뿐이다. 이와 같은 가정들을 고려하여, 일반적인 유체의 지배방정식인 Navier-Stokes 방정식으로부터 천수방정식을 유도할 수 있다.

점성을 고려하지 않을 경우, Navier-Stokes 방정식은 다음과 같이 표현된다:

∂ u → ∂ t = − u → · ∇ u → + f → ρ − ∇ p ρ, ∇ · u → = 0.

(1)

위 식에서 첫 번째 식은 유체의 속도에 대한 식이며, 두 번째 식은 비압축성에 대한 식이다. $u →$ 는 유체의 속도 필드, t는 시간, ρ는 밀도, $f →$ 는 외력, 그리고 p는 압력을 의미한다. 속도는 방향을 가지고 있으므로 3차원 벡터 (u,v,w)로 표현되며, 여기서 연직 방향의 성분 v에 대한 식만 살펴보면 다음과 같다:

∂ v ∂ t + u ∂ v ∂ x + v ∂ v ∂ y + w ∂ v ∂ z + 1 ρ ∂ p ∂ y = − g,

(2)

g는 중력가속도이며 연직방향의 성분에만 영향을 미친다. 위의 가정에 의해, 성분 v의 대류항들을 무시할 수 있으며, 이에 따라 중력가속도 g에 의한 압력의 관계식으로 간략화 된다:

1 ρ ∂ p ∂ y = − g .

(3)

y=h(h는 물 표면 지점의 높이)인 지점에서 p=0이므로 이를 p에 대한 식으로 다시 쓰면, 아래와 같다:

p (x, y, z, t) = p g (h (x, y, t) − y) .

(4)

수평 성분의 방정식에서 가정에 $∂ u ∂ y = ∂ w ∂ y = 0$ 의해 이며, 식(4)에 의해 압력 구배에 대한 항은 높이의 구배로 치환할 수 있다:

∂ u ∂ t + u ∂ u ∂ x + w ∂ u ∂ z + g ∂ h ∂ z = 0, ∂ w ∂ t + u ∂ w ∂ x + w ∂ w ∂ z + g ∂ h ∂ z = 0.

(5)

이처럼 연직 방향의 속도 변화를 무시하고 수평 성분의 속도로 관계식을 새로 정립함으로써 속도에 대한 3차원 미분방정식을 2차원으로 간략화 하였다.

일반적으로 유체의 표면은 레벨셋 값(level set value)[23] ϕ가 0인 지점으로 생각할 수 있으며, 천수의 경우 y=h인 지점에 해당한다. 따라서 ϕ(x,y,z)=y−h(x,z)로 치환이 가능하며, 이 값은 대류 방정식을 만족한다:

∂ ϕ ∂ t + u ∂ ϕ ∂ x + v ∂ ϕ ∂ y + w ∂ ϕ ∂ z = 0, − ∂ h ∂ t + u (− ∂ h ∂ x) + v (1) + w (− ∂ h ∂ z) = 0.

(6)

여기서, 수직 성분 v는 비압축성 조건을 이용하여 구할 수 있다. 비압축성 조건에 따라 y에 따른 v의 변화율은 아래와 같다:

∂ v ∂ y = − ∂ u ∂ x − ∂ w ∂ z .

(7)

v를 구하기 위해 필요한 유체와 지표면 사이의 경계조건은 $u → · n^= 0$ 이며, 지표면 b(x,z)에서의 법선벡터는 $v (− ∂ b ∂ x, 1, − ∂ b ∂ z)$ 에 비례하므로, 경계 조건에 적용하면,

v (x, y, z, t) = u ∂ b ∂ x + w ∂ b ∂ z − (∂ u ∂ x + ∂ w ∂ z) (y − b),

(8)

로 정리된다. 이제 식(8)을 식(6)에 대입한 후 정리하면 천수의 높이값에 대한 편미분방정식이 아래와 같이 유도된다:

∂ (h − b) ∂ t + u ∂ (h − b) ∂ x + w ∂ (h − b) ∂ z = − (h − b) (∂ u ∂ x + ∂ w ∂ z) .

(9)

여기서, 지면으로부터 수면까지의 거리(h-b)를 η로 치환하고 수평 방향의 속도, $v → = (u, w)$ 에 대한 식(5)과 함께 정리하면 최종적으로 천수방정식은 아래와 같다:

∂ v → ∂ t + (v → · ∇) v → = − g ∇ h . ∂ η ∂ t + (v → · ∇) η = − η ∇ · v → .

(10)

이 수식이 의미하는 바를 살펴보면, 각각의 물리량 η, $v →$ 는 수평 속도 $v →$ 에 따라 대류하고, 속도는 높이값 h의 구배에 의해 발생하며, η는 속도 필드의 발산(divergence)에 의해 변화하는 것을 알 수 있다. 이로써 풀어야 하는 공간 격자는 2차원으로 줄어들었다. 뿐만 아니라 압력을 미리 정해진 값으로 가정하고 있으므로 압력 분포를 알아내기 위한 선형 시스템을 풀 필요가 없는 것도 큰 이점이다.

3.2 천수방정식의 수치해법

천수방정식은 대류항(advection term)과 비대류항(non-advection term)으로 나뉘어지는데, 대류항은 Semi-lagrangian 기법으로 풀 수 있다.

식(10)의 좌변이 대류항인데 이를 다시 쓰면,

∂ v → ∂ t + (v → · ∇) v → = 0, ∂ η ∂ t + (v → · ∇) η = 0.

(11)

위 식은 속도와 물의 높이가 격자점 중심의 속도장을 따라 움직이는 것을 나타낸다. Semi-lagrangian 기법은 먼저, 속도장 위에 입자가 존재한다고 가정한다(Figure 1에서 격자 중앙에 빨간색 실선으로 표시된 원). 현재 격자의 중심에 입자들이 있고 속도장이 정의되어 있다면 그곳에 정의된 속도 벡터 값에 따라 그 입자들의 이전 타임 스텝에서의 위치를 역으로 추적하여 알아낼 수 있다(Figure 1에서 빨간 점선으로 표시된 벡터를 따라 이동한 지점의 원). 물론 역으로 계산된 입자의 위치는 격자의 중앙이 아닌 경우가 대부분이므로 인접한 네 개의 격자(X자 형태로 표시된 부분)를 기준으로 보간(interpolation)하는 과정이 필요하다. 이렇게 계산된 물리량을 현재의 값에 반영하여 대류항을 풀 수 있다. 대류 항을 푸는 알고리즘을 의사코드로 요약하면 다음과 같다.

Figure 1: Back-tracking and interpolation in the Semi-lagrangian method

물 표면 시뮬레이션을 위한 보존적 USCIP법

요약

Abstract

1. 서른

2. 관련 연구

3. 천수방정식의 수치해법

4. 천수방정식에서의 Constrained Interpolation Profile(CIP)법

5. 실험 결과

6. 결론 및 향후 연구

감사의 글

참고 문헌

<저자소개>