KR20050044147A - 연속체 유체역학 해석기법과 볼륨렌더링 기법을 이용한 실시간 유체유동 시뮬레이션 및 렌더링 방법 - Google Patents

Abstract

개시된 유체유동 시뮬레이션 및 렌더링 방법은, 연속체에 대한 유동방정식 해석방법에 기초하여 계산영역을 고정격자 형태의 유한요소로 분할하고 유체의 속도와 압력 등을 산출하여 자유표면의 형태를 산출하는 유체운동 시뮬레이션 단계와, 해당 유체를 통과하거나 유체의 표면에서 반사되는 빛의 성질을 이용하여 유체의 입체감을 표현하기 위한 시각적 변수들을 구하여 유체운동 시뮬레이션 단계에서 시뮬레이션된 유체에 대한 렌더링을 수행하는 볼륨 렌더링 단계를 포함한다. 이와 같은 방법에 의하면, 물리적 법칙을 엄격히 만족시키는 정확한 운동방정식을 사용하여 시뮬레이션을 수행하게 되므로, 기존의 유체효과 동영상 제작기법에 비해 높은 사실감의 영상물을 제작할 수 있으며, 작업자의 반복적인 수정작업이 불필요하여 전체 작업시간을 줄일 수 있고 복잡한 유동현상을 묘사하는 동영상도 실시간으로 구현할 수 있다.

Description

연속체 유체역학 해석기법과 볼륨렌더링 기법을 이용한 실시간 유체유동 시뮬레이션 및 렌더링 방법{A real-time simulation and rendering method for fluid flows using continuum fluid mechanics and volume rendering techniques}

본 발명은 유체유동 시뮬레이션 및 렌더링 방법에 관한 것으로서, 더 상세하게는 연속체 시뮬레이션 기법과 볼륨렌더링을 이용하여 유체의 복잡한 유동현상을 실시간으로 구현하는 방법에 관한 것이다.

컴퓨터그래픽 분야에 있어 유체의 자유운동을 표현하기 위한 방법의 하나로서, 비교적 단순한 역학적 방정식을 사용하여 유체의 질점 사이의 동역학적 상호작용 및 운동을 기술하는 입자 동역학 방법이 제안된 바 있다. 그러나 이러한 입자 동역학만으로 표현될 수 있는 범위가 상당히 제한적이며, 많은 사람들이 이런 방식으로 유체의 다양한 움직임을 표현하기에 많은 제약이 따른다는 사실에 공감하고 있다. 즉, 입자 방식에서는 기본적으로 전단응력(shear stress) 및 점성(viscosity), 회전성(vorticity) 등의 유체 고유의 다양한 성질을 효과적으로 표현하기 어려우며, 입자 사이의 충돌이나 벽면과의 충돌에 의해 유동경로가 급격히 변화할 경우 속도에너지를 상실한 입자들이 공간적으로 축적되어 총 부피의 보존법칙을 엄격히 만족시키지 못하는 등의 많은 문제점들이 존재한다.

유체 시뮬레이션에 있어 매우 중요한 또 한 가지 사항은 시뮬레이션 알고리즘의 안정성이다. 많은 경우에 연산에 적용되는 시간의 증분(increment)이 적절하지 못할 경우 연산알고리즘의 수학적 안정성에 따라 연산이 완료되지 못하고 부적절하게 중단되는 현상이 발생한다. 이러한 수학적 불안정성은 시뮬레이션 툴의 신뢰성을 떨어뜨리고 연산으로 얻어진 결과물의 시각적 개연성을 저하시키는 요인이 된다.

최근 컴퓨터 그래픽 분야에서도 입자 시뮬레이션의 표현의 제약 및 시각적 한계를 극복하기 위해 전통적인 유체역학분야의 기술을 활용하여 Navier-Stokes 방정식을 통한 연속체 시뮬레이션(Continuum Simulation) 방식으로 유체효과를 표현하는 연구들이 활발히 진행되고 있다. 1996년에 Nick Foster는 3차원 비압축성 Navier-Stokes 방정식을 적용하여 유동을 시뮬레이션 한 후에 가시화하였다. 경계 조건으로서 정지해있는 물체와 자유 표면을 표현하였고, 그 외에 부력 효과도 또한 포함하였다. 1999년 SIGGRAPH 에서, Stam은 어떠한 time step 값에 대해서도 Navier-Stokes 방정식을 비교적 정확하고 안정적으로 풀 수 있는 방법을 제시하여 이후 연구의 방향을 제시하였다. 이때 사용된 수치해석 기법은 기존의 입자 방식과는 달리 격자구조의 메쉬를 이용한 유한 차분법을 적용한 것이다. 1999년 Yngve 등은 역시 유한 차분법을 적용하여 충격파 문제와 폭발에 의한 화염 및 연무의 움직임을 해석하는 연구를 수행한 바 있다.

PDI studio의 Nick Foster 와 Stanford 대학의 Ronald Fedkiw 등에 의해 공동으로 진행되어 2001년 SIGGRAPH에 발표된 Practical Animation of Liquids의 연구에서는 입자 방식과 연속체 방식의 두 가지 방법을 상황에 맞게 혼용하고 level set 방식을 도입하였다. 이 방식은 유동해석에 적합한 물리적인 법칙을 적용하여 유체의 운동을 충실하게 재현하기 위해 시도된 것이다. 그러나 이들이 사용한 방법은 유체입자의 크기가 작은 많은 개수의 물방울이나 분수, 빗방울 등과 같이 수많은 미소입자로 구성된 유동문제를 해석하기 위해서 지나치게 많은 계산시간이 소요되며, level set 방식으로 표현된 유체의 자유표면 형상이 그들이 사용한 유동방정식의 수치적 확산에 따른 오류를 여과 없이 노출하게 되므로, 물리적 현상의 정확한 재현에 한계가 따르고 연산시간이 지나치게 오래 소요되어 실시간으로 신속한 계산이 요구되는 문제에는 적용하기에 어려움이 따른다.

렌더링 기법으로는 ray casting (광선 투영) 기법, blackbody radiation (흑체 복사) 기법 등이 적용되어 보다 사실적인 표현을 가능하게 하고 있다. 그 외의 기법으로는 물방울(blob) 개념의 유한한 크기를 가진 입자들을 사용한 입자 동역학을 활용하여 전체적인 움직임을 연산한 후 섬세한 텍스쳐는 매핑 기법으로 처리하는 등의 방식이 사용되고 있다.

사실적인 렌더링을 위해서 자연 현상의 렌더링에는 전역적 조명(global illumination) 방법이 많이 쓰였다. 그 중에서 Jensen에 의해 1998년에 제안된 photon mapping 을 이용한 방법은 최근의 유체 시뮬레이션에서 많이 쓰이고 있다. 이 방법을 통해 얻어지는 결과물의 시각적 품질은 매우 우수하나, 그에 반비례하여 수행 시간이 1 frame에 수분에서 수십 분이 걸리기 때문에 실시간 적용에는 적합하지 않다.

기체나 액체, 화염, 기후현상 등의 유체운동 시뮬레이션에는 기술적으로 다음과 같은 몇 가지 제약 요인들이 존재한다.

첫 번째 제약 요인은 유체운동의 역학적 시뮬레이션과 관련된 연산 시간이다. 전통적인 유체역학 분야에서는 정밀도 높은 유체 시뮬레이션을 위해서 Navier-Stokes 방정식을 이용한 계산 기법이 많이 사용되었다. Navier-Stokes 방정식을 기반으로 하는 유동해석 기술 중 자유표면의 운동을 모사하는 기술은 크게 고정격자법 (fixed grid method)과 이동격자법 (moving grid method)이 주로 사용되는데, 각각 장단점을 가지고 있다. 이동격자법은 움직이는 유체 경계면의 위치가 기하학적으로 명료하게 표현되는 대신 복잡한 형상의 영역을 움직이는 유체의 운동을 표현하기 위해서는 지나치게 많은 격자와 오랜 계산시간을 필요로 하는 단점이 있다. 이에 반해 고정격자법은 계산영역 전체를 초기에 미리 정해진 해상도로 분해하여 구성한 뒤 유체가 움직임에 따라 각 미소체적이 유체에 의해 점유되는 부피분율을 시간에 따라 변화시키는 방식으로 유체의 운동을 표현하게 된다. 이 방법과 더불어, 이동하는 대상물과 유체의 상호작용을 효과적으로 고려하기위해 고체분율(solid fraction)을 정의하여, 주어진 체적 내에서 고체가 차지하는 부피분율을 시간에 따라 변화시킴으로써 고체의 운동 및 그에 따른 유체의 교란현상을 보다 효율적으로 모사하는 방법이 개발되어 사용되고 있으나, 아직까지 실시간에 대응하기에는 연산속도가 느리다는 문제점을 가지고 있다.

두 번째 제약 요인은 렌더링 시간과 관련된 제약이다. 일반적인 포톤 사상(photon mapping)등을 이용한 전역적 조명(global illumination) 방법은 사실적인 렌더링을 위해 계산 시간이 상당히 많이 소요된다. 또한, 3D Max나 Maya 등의 상용 소프트웨어를 통한 시뮬레이션 데이터 렌더링도 역시 많은 시간이 소요된다. 이들 방식의 경우는 640×480 image의 5초 분량을 만드는데, 대략 20분 이상이 소요된다. 이러한 두 제약 요인을 해결하기 위해 여러 가지 기법들이 시도되고 있으며, 최근의 DirectX 기술에서는 파동전파(wave propagation), 텍스쳐 사상(texture mapping), 높이장(height-field) 등의 방식으로 실시간으로 유체효과를 표현하는 방법들이 개발되고 있다. 그러나 이러한 방식들은 임의의 자유표면을 표현하기에 부적절하여 평면적인 효과에 약간의 변화를 주는 수준에 머무르고 있다. Maya나 Max 등은 표준 툴의 플러그인 개념으로 Air와 Cloud 등과 같이 물리학을 기반으로 3차원의 대기와 구름 같은 환경효과를 만들어 주는 플러그인이 만들어져 있다. Air는 렌더링 시에 그늘진 하늘을 만들기 위해 라이트의 각도와 밝기를 사용하고 하늘의 색을 애니메이션하는 것을 담당한다. 또한 엷은 안개에서 폭풍우가 이는 안개는 물론 무지개 효과까지 다양한 시뮬레이션이 가능하다. Water는 사실적인 물의 표면을 생성한다. 이들 플러그인들에는 구름의 밀도, 뭉침 정도, 색 등의 파라미터들의 제어로 자연스러운 대기효과를 만들 수 있다. 이들 방식에서는 표현상의 제약은 많이 해소할 수 있지만, 사실적인 영상을 제작하기에는 데이터량이 상당히 방대하므로 실시간 랜더링에는 적합하지 않다. Nextlimmit, Inc.에서 개발한 RealFlow, RealWave, RealTracer는 스탠드어론 플러그인 형식으로 맥스, 마야, 라이트웨이브등에 소프트웨어에서 모두 사용이 가능한 유체관련 특수효과 프로그램이다. 최근 출시된 Maya 4.5의 fluid 모듈 역시 우수한 성능을 보여주고 있다. 그러나 이러한 프로그램들은 영상결과물의 품질이 우수하며 호환성이 뛰어나지만 기본적으로 particle 방식만을 채택하고 있으므로 모든 물리적 범위의 유동현상에 폭넓게 적용되기에 무리가 따르며, 실시간 렌더링을 지원하지 않는다. 실시간 가상현실 렌더링 툴(tool)로는 Java3D, OpenGL, DirectX 등의 다양한 표준들이 제작되고 있으며, 유체시뮬레이션을 적용한 일부 기능이 구현이 되고 있으나, 복잡한 유체유동이나 화염 등과 같이 부정형 유동을 실시간으로 효과적으로 시뮬레이션하기 위한 VR 솔루션은 아직까지 개발 단계에 머무르고 있어 그 품질은 매우 원시적인 수준이며 충분한 사실감을 구현할 수 없는 상태이다.

결과적으로 기존의 방법들은 시뮬레이션과 렌더링 시에 각각의 연산 시간 제약으로 인해 상호작용(interactive)을 할 수 없기 때문에 기존의 방법을 사용하여 다양한 효과를 만들어 내는 것이 시간이 많이 걸리고, 상호 작용을 한다고 하더라도 사실적인 영상을 얻어내기가 어렵다.

본 발명은 상기의 문제점을 감안하여 창출된 것으로서, 유체 복잡한 자유 운동을 실시간으로 신속하게 구현하고 물리적 법칙을 엄밀히 만족시켜 보다 사실감 높은 영상결과물을 제작할 수 있도록 개선된 유체유동 시뮬레이션 및 렌더링 방법을 제공하는데 그 목적이 있다.

상기의 목적을 달성하기 위한 본 발명의 유체유동 시뮬레이션 및 렌더링 방법은, 연속체에 대한 유동방정식 해석방법에 기초하여 계산영역을 고정격자 형태의 유한요소로 분할하고 유체의 속도와 압력 등을 산출하여 자유표면의 형태를 산출하는 유체운동 시뮬레이션 단계와, 해당 유체를 통과하거나 유체의 표면에서 반사되는 빛의 성질을 이용하여 유체의 입체감을 표현하기 위한 시각적 변수들을 구하여 상기 유체운동 시뮬레이션 단계에서 시뮬레이션된 유체에 대한 렌더링을 수행하는 볼륨 렌더링 단계를 포함하는 것을 특징으로 한다.

이하, 첨부된 도면을 참조하면서 본 발명에 따른 바람직한 실시예를 상세히 설명하기로 한다.

본 발명에서는 폭 넓은 물리적 영역의 유동현상 중 특히, 액체 및 기체, 화염 등의 불규칙한 운동을 일관된 방법으로 해석하고 이를 신속한 실시간 렌더링 환경에 접목하여 효율적으로 가시화한 실시간 시뮬레이션 프레임워크(framework)를 제안하고자 한다. 본 방법은 인터렉티브에 대응할 수 있는 속도와 일반적인 실시간 응용분야에 활용하기에 적절한 품질의 이미지를 보여주며, 특히 유체의 급격한 속도변화나 국소적인 정체현상 등 극단적인 경우에 있어서도 물리적 법칙에 입각한사실적인 영상물을 안정적으로 도출하도록 해준다.

본 발명의 유체 동영상 구현방법은 도 1에 나타낸 바와 같이 크게 3단계로 나뉘어 진행되는데, 연속체 시뮬레이션 단계(S1)와, 그 시뮬레이션 된 연속체의 입체감을 나타내기 위한 데이터를 산출하여 볼륨렌더링을 수행하는 단계(S2) 및, 보다 사실감 있는 시각적 효과를 부여하기 위해 사용자로 하여금 각종 파라미터를 조정할 수 있도록 해주는 유동효과 조정단계(S3)가 그것이다.

우선, 시뮬레이션 단계(S1)에서는 유한 요소에 의해 유동을 표현하는 유체 역학을 기반으로 시뮬레이션 하였다. 이 방식은 입자 동역학 방식의 단점인 제한된 개수의 입자 혹은 격자점의 데이터의 불연속성으로 인한 사실감의 결여 없이 부드럽고 자연스럽게 연결되는 표면을 표현할 수 있다. 또한 연속체 시뮬레이션에 사용되는 유한요소법 기반의 유체해석용 격자 구조와 렌더링에 사용되는 voxel 개념의 3차원 격자구조는 그 기하학적 구조가 매우 유사하여 시뮬레이션을 통해 얻은 데이터를 손쉽게 렌더링에 이용할 수 있었다.

그러나 연속체 유동방정식을 이용해 표현하고자 하는 유체는 일반적으로 응집된 덩어리의 형태를 가지고 있으므로 한 연속체 내에 많은 개수의 단위격자가 포함된다. 이러한 규칙적 격자구조를 3차원적으로 렌더링할 경우 관찰자의 시점에 따라 각 층마다 특유의 반복패턴이 나타나 시각적 사실감을 저하시키는 문제점이 발생한다. 이러한 단점을 극복하고자 제2단계인 볼륨렌더링 단계(S2)에서는 복셀(voxel) 구조에 기초한 3차원 실시간 예비 적분 볼륨 렌더링(pre-integrated volume rendering) 방법을 사용하였다. 이와 같이 본 발명에서 응용한 볼륨 렌더링 방법은 Levoy에 의해 1990년에 제안된 ray casting 방법으로부터 시작되었다. 이 방법은 많은 양의 볼륨 데이터를 계산하기 때문에 속도가 느리지만 공간 도약법과 조기 광선 중단법을 통해 가속화를 할 수 있다. 최근의 그래픽스 하드웨어의 발달로 실시간 볼륨 렌더링을 위한 하드웨어 텍스쳐를 이용한 가속법이 많이 연구되고 있다. Engel은 2001년에 하드웨어 텍스쳐를 이용한 가속화 방법에 더하여 OTF의 문제점을 보완한 상기의 예비 적분 (pre-integration) 방법을 제안하였다. 본 발명의 렌더링 방법은 기본적으로 이 Engel의 방법에 기초하고 있다.

그리고, 유동효과 조정단계(S3)에서 각종 파라미터의 조정과 다중 OTF(Opacity Transfer Function, 불투명도 전이 함수)의 사용을 통해서 다양한 종류의 유동을 거의 실시간에 표현할 수 있게 된다.

이하, 상기한 시뮬레이션 단계(S1)부터 차례로 설명하기로 한다.

본 발명에서는 상기한 바와 같이 유한 요소법으로 유동을 표현하는 유체 역학을 기반으로 하여 시뮬레이션을 수행한다. 일단, 자유 표면을 갖는 유체는 경계 부분이 미리 정해지지 않고, 해석과정에 따라 결정되므로 경계가 움직이는 문제점이 있다. Navier-Stokes Equation과 마찬가지로 이 문제는 매우 비선형적인 특성을 갖는다. 영역의 연속적인 변화와 경계면의 물리적 비연속성은 이 문제를 다루기 매우 어렵게 한다.

이와 같이 움직이는 경계면 문제를 해결하기 위한 방법으로는, Lagrangian방법 (도 2a 참조), Eulerian방법(도 2b 참조), 그리고 두 방법을 결합한 ALE(Arbitrary Lagrangian Eulerian) 방법이 있다. Lagrangian 방법에서는 좌표계가 유체와 같은 속도로 움직이고 자유표면은 메쉬와 일치한다. 그 결과, 경계 면이 정확하게 얻어진다. 그러나, 왜곡이 심하고 형상이 복잡한 유체는 정확한 경계 면을 얻기가 어렵다. 이런 경우, 계산 시간이 많이 소요되지만 왜곡된 메쉬를 재접합하는 작업이 불가피하다.

Eulerian방법은 좌표계는 고정되어 있고 유체가 좌표계 내에서 움직인다. 또, 처음에 만들어진 메쉬가 전체 계산시간동안 사용된다. 이 방법을 사용하면, 움직임이 큰 유체도 비교적 정확하고 효과적으로 다룰 수 있다. 또한, 여러 경계면을 다룰 수 있고 2차원 코드를 3차원으로 확장하기도 쉽다. 반면에, 경계면이 격자점과 일치하지 않고 격자점들 사이에서 움직이므로 정확하게 경계면을 얻어낼 수 있는 특별한 방법이 필요하다. Eulerian 방법은 고정격자법(fixed grid method), 이동격자법(moving grid method), 매핑 방법(mapping method)의 3가지 유형으로 나뉜다.

고정격자법은 전체 연산영역을 미리 적당한 크기의 격자구조로 분할한 후 각 격자내에 유체의 질량이 유입될 때 유입된 유체의 부피와 격자의 고유부피 사이의 비율을 이용하여 유체의 공간적 점유율을 스칼라함수로 정의하고, 그 값을 시간에 따라 변화시켜 유체가 점유한 영역의 전체적인 형상을 추적하는 방법이다. 대표적인 기법으로는 VOF(Volume of Fluid) method(도 3a), MAC(Marker And Cell) method(도 3b), SLIC(Simple Line Interface Calculation) method, 그리고 Young의 방법 등이 있다.

이동격자법은 경계면과 일치하도록 메쉬를 조정함으로써 경계면을 정확하게 위치시키는 방법이다. 계산이 진행됨에 따라 자유표면에 위치한 유체 질점의 속도가 변화하게 되며, 이 데이터를 사용하여 유동표면의 전진에 따른 연산영역의 변화된 형상을 기하학적으로 계산하여 매 시간마다 새로이 연산영역을 정의해 주는 방법이다.

매핑 방법은 물리적인 영역이 매핑 함수에 의해 간단한 계산 영역으로 변형된다. 이 방법은 경계면을 정확하게 풀 수 있으나, 전체 형상을 적절한 수학적 함수의 형태로 표현할 수 있어야 하므로 복잡한 형상에 대해 이 방법을 적용하기는 매우 어렵다. 따라서 복잡한 경계면을 가진 유동에 대해서는 적절한 매핑 함수를 찾아내기가 어렵다.

Lagrangian과 Eulerain방법을 결합한 Arbitrary Lagrangian-Eulerain(ALE)방법은 위의 두 방법의 단점은 피하고 장점만 결합하였다. 이 방법에서 메쉬의 속도는 유체의 속도와 독립적이나 경계면의 메쉬는 경계면을 쉽게 위치시키도록 유체와 같은 속도로 움직인다.

본 발명에서는 상기 고정격자법 중 VOF 기반의 방법을 사용하여 불규칙한 유동을 추적하였고, 자유표면을 재구성하기 위해서는 baby-cell 방법을 이용하였다. 이 방식은 김(Kim)에 의해 제안되었으며, 수학식 14 이후에 설명하기로 한다.

이하, 비압축성 유체의 해석 알고리즘을 설명하기로 한다.

먼저, 지배방정식을 이용하는 방법에 대해 설명한다. 비압축성 유체에 대한 지배 방정식은 수학식 1과 같은 연속체 방정식과 수학식 2의 Navier-Stokes Equation이 있다.

여기서, 는 속도 벡터, p는 압력, 는 중력 또는 부력, 는 viscous stress tensor, 는 밀도이다.

그리고, 구성방정식은 다음과 같이 표현된다.

는 연신율 텐서(strain rate tensor)이고 는 전체 응력 텐서(stress tensor)이다. 는 점성이다.

속도를 위한 경계 조건(boundary condition)은 의 필수경계조건(essential boundary condition)과 의 정지마찰의 경계조건(traction boundary condition)으로 구성된다. 수학식 5는 경계 조건을 표현한 식이다.

상기 수학식 5의 두 조건은 상호배제적이다. 는 의 속도 벡터이고 는 의 트랙션(traction) 벡터이다. 는 의 단위 normal 벡터이다.

전체 영역 에 대한 초기 조건은 다음과 같다. 에서 와 에서 의 조건을 만족하면 가 해답이 존재하기 위한 초기 조건이 된다. 만약 가 null이라면, 속도를 위한 경계 조건은 반드시 다음의 수학식 6 조건을 만족해야 한다.

여기서 는 속도의 초기 값이다.

유한요소법(Finite Element Method)은 복잡한 형상의 전체 영역을 유한개의 간단한 유한 요소로 나누어 단순한 기하학적 요소의 집합으로 다루며, 지배 방정식에 대한 경계 조건도 다루기 쉽도록 변형시켜 해석하는 수치해석 기법이다. FEM은 공학과 수학적 물리학 등 많은 영역에서 널리 사용된다. 본 발명에서는 FEM을 사용하여 지배 방정식을 해석하였다.

비압축성 Navier-Stokes 방정식을 위한 FEM 중에는 mixed method, penalty method, fractional step method가 있다. mixed method에서 속도와 압력은 동시에 해석된다. 즉 질량과 운동량 보존이 동시에 만족된다. penalty method를 사용할 경우 연속체 방정식과 운동량 방정식 사이에 적절한 수치적 가정을 도입하여 양쪽을 모두 적당한 오차한계의 범위 이내로 수렴시키는 과정이 개입된다. 이 방법을 사용할 경우 전체의 지배방정식이 속도만에 대한 방정식으로 간략화되며, 압력은 필요한 경우에 추가적으로 계산된다. fractional step method는 운동량 보존 식과 질량 보존 식을 순차적으로 모두 만족시킨다.

다음으로 Explicit element-by-element fractional step method라는 방법이 있다.

이 방법은 나카야마(Nakayama)와 모리(Mori)에 의해 제안되었으며, 김(Kim) 등에 의해 응용되었다. 반 이산화 운동량 방정식(semi-discretized momentum equations)은 다음과 같이 표현될 수 있다.

그리고, 연속체 방정식은 다음과 같이 이산화될 수 있다.

여기서, u와 p는 각각 node의 속도와 압력 벡터이다. M, C(u), K, H는 질량, 대류, 발산, 압력 변화도 행렬이다. F는 외부의 힘과 경계 조건을 포함하는 힘의 벡터이다. 는 발산 행렬을 의미하며 N은 요소의 개수이다. Euler 방법이 반이산화된 운동량 방정식에 다음과 같이 적용될 수 있다.

는 매개 속도를 나타낸다. 상기 수학식 9에서 삼선형 모양 함수가 속도 변수에 적용되고 단위요소 모양 함수는 압력 변수에 적용된다. 매개 속도 는 발산하지 않는다고 볼 수는 없으나 압력 을 가지고 질량 보존 법칙을 만족하도록 수정될 수 있다. 압력은 요소의 발산을 없애도록 반복적으로 요소요소를 일소함으로서 조정된다. 속도는 압력에 의해 수정된다. 속도와 압력 수정 단계는 다음과 같이 요약될 수 있다.

처음에 와 의 값은 각각 매개 속도 와 이다. 단위 볼륨 당 발산은 다음과 같이 계산된다.

n+1과 (k)는 time step과 반복 회수를 말한다. 는 요소의 볼륨을 말하고, 는 전체 압력 변화량 행렬 H에 상응하는 요소의 값이다. 만약 의 크기가 수렴 범위 이내라면, 유체는 지역적으로 비압축성이고 더 이상 조정이 필요하지 않다. 만약 을 넘는다면 압력은 발산하지 않는다는 조건을 만족하도록 조정될 수 있다. 속도는 수학식 11과 같이 발산하지 않도록 압력에 의해 수정될 수 있다.

이 방법은 모멘텀 보존 식과 속도 수정 단계 둘 다에서 row-sum mass lumping 방법을 적용하였다. 본 발명에서 은 5×10^-4에서 1×10^-5까지다. 이 방법은 수렴을 위해 여러 번 반복연산이 필요하지만 전체적인 성능이 우수하다. 움직이는 자유 표면을 가진 유체는 일반적으로 복잡하고 매우 계산시간이 많이 소요된다. explicit element-by-element step방법은 비압축성 유체를 계산하는데 매우 효율적인 알고리즘이다. 본 발명에서는 효율성을 위해 이 방법을 채택하였다.

다음으로, 자유 표면을 가진 유체의 해석 알고리즘에 대해 설명하기로 한다.

일반적으로 자유 표면은 함수가 매우 복잡하고 유체가 매우 빠르게 움직이는 특성을 가진다. 본 발명에서는 복잡한 유체의 움직임을 위해 상기한 바와 같이 Lagrangian, Eulerian, ALE 방법에 기반 한 많은 방법들 중에서 VOF 방법을 이용하였다. 그리고, 본 발명에서는 다음을 가정하였다. 첫째, 시스템은 두 개의 혼합되지 않는 유체로 구성된다. 둘째, 혼합되지 않는 두 유체의 물리적 성질은 내부적으로 균질하다. 셋째, 유체는 비압축성이고 점성이 있으며 관성적 영역은 난류현상이 발생하기 직전까지의 범위를 포함한다. 넷째, 표면 장력의 효과는 무시된다.

VOF방법은 Hirt와 Nichols에 의해 처음 개발된 것으로 움직이는 자유 표면을 가진 유체를 해석하기 위한 고정좌표 시스템에 기반한 대표적 알고리즘이다. 유체 f(x, y, z, t)의 볼륨은 수학식 12와 같이 한 단위요소 안에 유체의 비율로 정의된다.

한 단위요소가 완전히 유체로 채워지면, 유체의 볼륨은 f=1이고 그 단위요소는 주된 유체의 영역으로 간주한다. 단위요소가 비어있으면( f = 0 ), 그것은 빈 영역에 속하고 유체의 계산에서 배제된다. 단위요소는 유체의 볼륨이 0에서 1사이일 때 즉, 0 < f < 1일 때 자유 표면 위에 있는 것으로 생각한다. 도 4에 보여지듯이, 자유 표면의 계산 영역은 변수 f를 사용하여 쉽게 표현될 수 있다.

경계 면의 움직임은 f에 대한 다음의 수학식 13으로 표현된다.

이 식은 오직 비정상항과 대류항만으로 이루어져 있다. 그러므로 이 식을 직접 이산화하여 수치적 연산을 수행할 경우 해석은 유체 표면의 경계가 흐려지면서 심각한 발산을 유발할 수 있다. VOF 방법은 복잡한 움직임의 유체유동을 효율적으로 다룰 수 있으나, 경우에 따라서는 이러한 수치적 발산을 어떻게 줄일 것인가 하는 문제가 발생한다. 안정된 수치해석 결과를 얻기 위해서는 수치적 발산을 피할 수 있어야 하며 동시에 유체 표면의 경계가 흐려지지 않고 이동하여야 한다. 발산을 억제할 수 있도록 고안된 알고리즘으로는 SLIC (Simple Line Interface Calculation) method(도 5a 참조), SOLA-VOF method(도 5b 참조), Young's method(도 5c 참조), FLAIR(Flux Line-Segment Model for Advection and Interface Reconstruction) method(도 5d 참조), PLIC (Piecewise-Linear Interface Calculation) method(도 5e 참조)등이 있다. 이들은 유체 표면의 형상을 적절히 분할하여 재구성하여 각 유한요소 단위로 상황에 맞게 적절히 자유표면을 움직이는 알고리즘으로서, 계산의 정확도를 높이고 발산을 막는다.

자유표면을 전진시키는 기법에 있어 처음 고려할 문제는 자유 표면의 기하학적 재구성과 관련된 wet-out fraction을 어떻게 찾을 것인가 이다. wet-out fraction 과 을 계산하는 알고리즘은 김(Kim)의 baby-cell 방법을 기초로 한다. 유체 표면의 형상을 재구성하는 알고리즘 중에는 상기한 PLIC-type(도 4e)이 있는데 이 방법은 경계면을 normal vector가 f의 기울기(gradient)인 선이나 면으로 가정한다. baby-cell 방법은 PLIC-type을 사용하여 자유 표면을 표현한다. 이 방법에서 방향벡터 r은 수학식 14와 같이 정의된다.

와 는 각각 i 번째 요소의 j 번째 면의 유체의 볼륨과 단위요소의 볼륨을 의미한다. 또한, 는 요소의 경계 에 대한 단위 normal vector이고 는 i 번째 요소의 전체 영역 중에서 j 번째 면의 비율이다. 수학식 3과 9에서 볼 수 있는 바와 같이, 방향 벡터는 이웃한 유한요소 내의 유체의 체적평균을 구하여 결정된다. 그리고, 도 6에서 보듯이 방향벡터는 유체의 부피가 큰 방향을 나타내는 단위 벡터이다. 방향 벡터 r은 자유 표면의 normal 벡터 대신 사용될 수 있다. 이 벡터는 대수적인(algebraic) 방법으로 정의되므로 비구조적 메쉬(unstructured mesh)에서도 이 벡터를 계산하는 일은 매우 용이하다. 방향 벡터는 도 7에 나타낸 바와 같이 다음의 방법으로 wet-out fraction을 결정하는데 사용된다. 첫째, 방향 벡터 r에 수직되고 유한요소의 중심점을 통과하는 평면을 만든다. 둘째, 각 유한요소는 같은 체적을 가지는 baby-cell들로 분할된다. 셋째, baby-cell은 방향벡터 r의 양의 방향의 평면에서 가장 먼 cell부터 시작해서 채워진다. 그리고 나서, wet-out fraction 가 계산될 수 있다. FEM에서 좌표계 변형 속성을 사용하여 위의 계산이 수행된다. 이러한 시뮬레이션 과정을 통해 유체의 유동이 해석되고 각 셀에서의 속도와 압력 등이 구해지게 된다.

다음으로는, 볼륨렌더링을 수행하는 단계(S2)로 들어간다. 볼륨 렌더링은 유체 입자가 빛을 투과할 수 있는 물체로 가정하고 밝기, 불투명도, 색깔 등을 산출하여 입체감을 부여하는 과정이며, 먼저 상기 연속체 시뮬레이션으로부터 볼륨 데이터를 산출한다. 시뮬레이션 결과를 볼륨 렌더링에 이용하기 위해서는 시뮬레이션의 격자를 도 8에 도시된 바와 같이 볼륨 렌더링의 격자로 매핑하는 과정이 필요하다. 움직이는 유체의 해석 알고리즘인 VOF 방법에서 한 단위요소의 볼륨에 대한 유체의 볼륨의 비율인 f값은 단위 요소에 유체가 얼마나 채워져 있는가에 대한 값이다. 이 값을 도 8과 같이 단위 체적을 대표하는 복셀(voxel)의 데이터로 매핑하면 볼륨 렌더링을 위한 데이터가 만들어지게 된다. 여기서, 복셀은 볼륨의 단위가 되는 각 셀의 점을 의미한다.

연속체 유한요소 해석에 사용되는 격자 구조와 볼륨렌더링의 복셀 구조는 공통적으로 육면체의 배열을 가지므로 그 형태가 매우 유사하다. 유한요소격자와 복셀의 셀 간격이 동일할 때는 각 유한요소의 fill fraction (수학식 12에서 정의된 f) 값을 그대로 복셀 데이터에 일대일 대응시킬 수 있다. 이때 각 유한요소의 절점에 그 점을 둘러싼 미소체적의 질량을 집중질량으로 취급하여 복셀 구조의 각 절점에 데이터를 부여한 후 이를 그대로 렌더링에 사용할 수 있다.

유한요소의 격자구조와 복셀데이터의 격자구조가 일치하지 않을 경우 각 유한요소 절점의 위치에 할당된 집중질량은 복셀데이터 격자 내에서 특정 육면체 내부의 한 위치를 차지하게 된다. 이때 유한요소절점의 집중질량이 차지하는 유효부피가 복셀데이터 격자내의 한 육면체 내에서 얼마만큼의 부피를 점유하는지를 계산하여 이 부피분율을 그 복셀 절점의 부피분율로 사용하게 된다. 경우에 따라 집중질량의 유효부피에 상당하는 반경이 한 개의 복셀 크기를 초과할 경우, 그 잉여 부피는 인접 복셀에 전가되어 더해지게 된다. 이러한 방식으로 전체의 유한요소 절점에 대해 유체가 차지하는 fill fraction을 복셀 데이터에 일대일 대응시키고 이어서 렌더링 과정을 거쳐 유체의 분포를 가시화하게 된다.

기본적인 볼륨 렌더링 식은 예비 적분(pre-integration)을 이용한 다음과 같은 형태로 주어진다.

I: 최종밝기값

λ: 시점으로부터의 거리,

x: 3차원 공간 좌표,

s: 연속적인 스칼라 장,

: 기본 칼라,

: 감쇄도,

D: 시점으로부터의 최대 거리.

여기서, d = D/n이라고 할 때, 수학식 15의 지수 함수 내부는 다음의 수학식 16과 같이 근사될 수 있다.

여기서 는 i번째 광선 부분의 불투명도를 나타낸다.

결국, i 번째 광선 부분(ray segment) 의 불투명도 는 다음의 수학식 17과 같이 근사될 수 있다.

그리고, 수학식 15에서 i번째 광선 부분에서 방출되는 색깔 는 수학식 18로 근사될 수 있다.

결국, 수학식 15의 적분은 다음의 수학식 19와 같이 정리될 수 있다.

즉, i번째 광선 부분의 최종밝기값은 0부터 i-1까지의 불투명도를 곱한 값과 i번째 광선 부분의 색상값을 곱해서 구해진다. 그리고 0부터 n번째 까지의 광선 부분의 최종밝기값을 모두 합산한 값이 I가 된다.

따라서, 후전방 합성 알고리즘(back-to-front compositing algorithm) 은 다음의 수학식20과 같이 계산된다.

은 n번째 광선 부분에서부터 i번째 광선 부분까지의 밝기의 총합을 나타낸다.

또한, 는 로 교체될 수 있으므로, i번째 광선 부분에서 방출되는 색깔 는 수학식 21과 같이 근사될 수 있다.

위 식을 이용하면, 수학식 15의 적분은 다음의 수학식 22과 같이 더 일반화된 형태로 표시될 수 있다.

즉, 수학식 19의 는 색깔과 불투명도가 곱해진 값이므로 풀어서 로 쓴 것이다.

i번째 부분(segment)에 대해서 segment의 앞쪽의 스칼라 값 (밀도 값) s_f:=s(x(id))가 되고, segment의 뒤쪽의 스칼라 값 s_b:=s(x((i+1)d)) 가 된다. 이 값을 이용하면, i번째 segment의 불투명도는 수학식 23과 같이 근사 될 수 있다.

w는 적분 파라미터를 의미한다.

결국 는 s_f, s_b, d의 함수가 된다. 따라서, i번째 광선 부분에서 방출되는 색깔 는 다음의 수학식 24와 같이 계산될 수 있다.

여기서 를 정의하면, 를 수학식 25과 같이 T를 사용해서 나타낼 수 있다.

비슷한 방법으로 를 정의하면, 를 K를 사용해서 수학식 26처럼 나타낼 수 있다.

이와 같은 수식들을 이용하여 최종 밝기와 색깔 및 불투명도 등 볼륨을 표현하기 위한 데이터를 산출되며, 이 데이터를 바탕으로 실시간으로 렌더링이 진행되는 것이다. 도 9a 내지 도 9c는 소정 용기 안에 유체가 부어지는 장면을 묘사한 렌더링 결과를 보인 것이다.

다음으로, 사용자로 하여금 파라미터를 조정할 수 있도록 입력 환경을 제공하고, 그 입력값을 반영하여 렌더링 결과물에 다양한 유동 효과를 부여하는 유동효과 조정단계(S3)로 들어가게 된다. 이때 사용되는 시뮬레이션 파라미터에는 유체의 밀도, 공기의 점도, 물의 점도 등이 있고, 렌더링 파라미터로는 하기의 수학식 18에 표시된 퐁(phong) 조명 모델식의 주변조명상수 k_a, 분산조명상수 k_d, 반사조명상수 k_s, 지수항 n 등이 있다. 예를 들어 n값을 조정하면 반사 조명의 밝기 값을 변화시켜서 반짝거림의 정도를 표현할 수 있다.

여기서, : 물체의 수직 벡터 (도 10a 참조), : 빛의 방향 벡터, : 빛의 반사벡터, : 시점과 물체가 이루는 방향벡터, I : 밝기, n : 지수항, (도 10b 참조)를 각각 의미한다.

따라서, 이와 같은 파라미터를 사용자가 조정하면 그에 상응하는 시뮬레이션 결과가 실시간으로 얻어지게 된다.

이와 같은 본 발명의 유체유동 시뮬레이션 및 렌더링 방법은, 상기와 같이 실시간으로 다양한 파라미터를 조정하면서 유체의 세밀한 이미지 품질을 신속하게 얻을 수 있으므로, 다양한 가상현실 가시화 기술에 유효하게 사용될 수 있다.

예를 들면, 유체나 화염 등의 복잡한 자연현상에 대한 가상현실(Virtual Reality) 가시화 기술은 3차원 입체영상(Stereoscope) 컨텐츠, 디지털 문화재 복원, 가상현실 여행 및 관광, 자연 및 체험, 교육 및 학습 등의 다양한 분야에 효과적으로 적용될 수 있다. 또한 이러한 가상현실 가시화 기술은 문화 컨텐츠 분야뿐만 아니라 교육용, 의료용, 공학용 등의 전문 분야에도 효과적으로 적용될 수 있다. 폭우나 짙은 안개 속의 항공기의 조종 연습 등과 같은 특수 교육 목적의 가상현실 컨텐츠에도 효과적으로 적용될 수 있다. 특히 의료 및 공학용 데이터들은 그 자료가 방대하고 정밀도가 요구되는 반면, 그 데이터를 사실적으로 가시화하기 위한 시각화 기술이 미흡한 실정이다. 각종 가전 및 전자제품을 비롯하여 고속 서버급 컴퓨터의 냉각, 순환을 위해 사용되는 냉각기술과 유동제어 기술, 자동차 및 항공기 등의 몸체 외부의 공기 유동이나 선박 주위의 유체유동 등의 해석은 성능개선을 위한 최적화 설계에 핵심적인 항목이다. 이러한 공학적 해석분야에 적용하여 그 결과를 효과적으로 시각화하기 위해서는 가상현실기법을 도입한 사실적인 시각화 기술이 필수적이다. 이외에도 체내 각 기관 내부의 체액의 유동이나 각 기관의 섬모운동 가시화, 폐 속의 공기나 담배연기의 흡입, 혈액이나 소변의 흐름, 소화기능 등의 가시화를 비롯한 다양한 의료, 생명공학 분야를 비롯하여 여러 가지 분야에 폭넓게 사용될 수 있다. 이 외에도 web 3D 렌더링이나 온라인 방송, HDTV (high-definition television)용 인터액티브 컨텐츠 등의 분야에서도 다양한 활용 가능성을 예측할 수 있다.

본 발명의 효과를 확인해보기 위한 예로서, 본 발명의 방법과 기존의 소프트웨어인 Nextlimit사의 realflow를 사용하여 용기에 물이 부어지는 약 10초 분량의 장면을 구현해보았다. 그 결과 비슷한 정밀도의 동영상 결과물을 도출하는데 있어 realflow를 사용한 경우에는 렌더링까지 끝난 결과물이 얻기까지 여러 가지 반복작업과 연산 및 렌더링을 위해 총 수십 시간의 시간이 소요되는데 비해, 본 발명의 방법은 필요한 파라미터만 정확히 입력되면 거의 실시간으로 수 십초 이내에 결과물을 얻을 수 있었다. 특히 유동의 형상이 복잡한 경우 기존의 소프트웨어로 작업하면 유동의 방향이 크게 선회하거나 급격한 정체현상이 발생하는 부분에서 물리적 제한조건이 만족되지 않아 작업자가 일일이 이를 수작업으로 보정하여야 하고 만족스러운 결과물의 제작이 불가능한 경우도 발생하였다. 그러나 본 발명의 방식을 사용할 경우 적절한 초기조건과 경계조건만을 입력하면 그 이후의 과정은 기하학적 형상과 유체의 물성, 외력 등의 요소에 의해 엄격한 물리적 법칙에 의해 자동적으로 모든 유체의 움직임이 계산되므로 동영상제작에 필요한 시간이 획기적으로 단축될 수 있다.

상술한 바와 같이 본 발명의 유체유동 시뮬레이션 및 렌더링 방법은, 다음과 같은 효과를 제공한다.

첫째, 연속체 시뮬레이션과 볼륨 렌더링을 채용함으로써 물리적 법칙에 더욱 충실한 동영상을 구현할 수 있으며, 수치해석 기법을 개선하여 연산시간을 줄임으로써 실시간으로 유체 동영상을 구현할 수 있다.

둘째, 사용자의 파라미터 입력을 받아서 다양한 동영상 이미지를 즉시 변화시킬 수 있기 때문에, 사용자가 원하는 동영상을 보다 신속하고 정확하게 구현할 수 있도록 해준다.

본 발명은 상기에 설명되고 도면에 예시된 것에 의해 한정되는 것은 아니며 다음에 기재되는 청구의 범위 내에서 더 많은 변형 및 변용예가 가능한 것임은 물론이다.

도 1은 본 발명에 따른 유체유동 시뮬레이션 및 렌더링 방법을 나타낸 플로우 챠트,

도 2a 내지 도 3b는 움직이는 유체 경계면을 해석하는 방법의 예를 보인 도면,

도 4는 자유 표면을 가진 유체를 해석하기 위한 고정좌표계의 예를 보인 도면,

도 5a 내지 도 5e는 유체의 표면 형상을 재구성하는 알고리즘의 예를 보인 도면,

도 6은 자유 표면의 각 요소에서의 방향 벡터를 보인 도면,

도 7은 도 6의 방향 벡터를 사용하여 wet-out fraction을 구하는 과정을 보인 도면,

도 8은 시뮬레이션의 격자를 볼륨 렌더링의 격자로 매핑하는 과정을 보인 도면,

도 9a 내지 도 9c는 본 발명의 유체유동 시뮬레이션 및 렌더링 방법에 의해 소정 용기에 부어지는 유체가 표현된 것을 보인 도면,

도 10a와 도 10b는 퐁 조명 모델의 각 변수를 설명하기 위한 도면.

Claims (4)

1. 연속체에 대한 유동방정식 해석방법에 기초하여 계산영역을 고정격자 형태의 유한요소로 분할하고 유체의 속도와 압력 등을 산출하여 자유표면의 형태를 산출하는 유체운동 시뮬레이션 단계와,

   해당 유체를 통과하거나 유체의 표면에서 반사되는 빛의 성질을 이용하여 유체의 입체감을 표현하기 위한 시각적 변수들을 구하여 상기 유체운동 시뮬레이션 단계에서 시뮬레이션된 유체에 대한 렌더링을 수행하는 볼륨 렌더링 단계를 포함하는 것을 특징으로 하는 유체유동 시뮬레이션 및 렌더링 방법.
2. 제1항에 있어서,

   상기 볼륨 렌더링 단계의 입체감을 표현하기 위한 변수들은 유체의 밝기, 색깔, 불투명도 등이며, 이들 변수들은 예비 적분(pre-integration) 방식을 이용하여 산출되는 것을 특징으로 하는 유체유동 시뮬레이션 및 렌더링 방법.
3. 제1항에 있어서,

   상기 유체유동 시뮬레이션과 볼륨 렌더링 단계의 연산을 위한 파라미터들을 사용자가 입력할 수 있도록 제공하고, 그 입력된 값을 실시간으로 유체 동영상에 반영하는 유동효과 조정단계를 더 포함하는 것을 특징으로 하는 유체유동 시뮬레이션 및 렌더링 방법.
4. 제3항에 있어서,

   상기 파라미터들에는 유체의 밀도와 점도 및 조명 상수 중 적어도 하나가 포함되는 것을 특징으로 하는 유체운동 시뮬레이션 및 렌더링 방법.