자유표면 유동해석을 위한 WMLS 기반 입자법 기술 개발

Development of WMLS-based Particle Simulation Method for Solving Free-Surface Flow

Article information

J. Ocean Eng. Technol. 2014;28(2):93-101
Corresponding author Jong-Chun Park: +82-051-510-2480, jcpark@pnu.edu
Received 2014 February 03; Accepted 2014 April 15.

Trans Abstract

In general, particle simulation methods such as the MPS(Moving Particle Simulation) or SPH(Smoothed Particle Hydrodynamics) methods have some serious drawbacks for pressure solutions. The pressure field shows spurious high fluctuations both temporally and spatially. It is well known that pressure fluctuation primarily occurs because of the numerical approximation of the partial differential operators. The MPS and SPH methods employ a pre-defined kernel function in the approximation of the gradient and Laplacian operators. Because this kernel function is constructed artificially, an accurate solution cannot be guaranteed, especially when the distribution of particles is irregular. In this paper, we propose a particle simulation method based on the moving least-square technique for solving the partial differential operators using a Taylor-series expansion. The developed method was applied to the hydro-static pressure and dam-broken problems to validate it.

1. 서 론

자유표면의 대 변형이나 동적 거동에 대한 해석은 일반적으로 격심한 비선형 자유표면 유동을 동반하기 때문에 해양공학 분야의 유동 시뮬레이션에 있어서 가장 다루기 어려운 문제 중 하나이다. 이러한 문제들을 다루기 위하여 기존에 주로 사용되는 방법으로는 오일러(Eulerian) 접근법인 격자법이 있다(Hirt and Nichols, 1981; Takewaki and Yabe, 1987; Sussman et al., 1994; Miyata and Park, 1995). 하지만 이러한 격자법의 경우에는 자유표면 근처에서 격자의 정도 높은 조밀성이 요구되며, 대류항에 의한 수치적 확산에 의해 질량보존이 완벽하게 만족되지 않는다. 한편, 라그란지(Lagrangian) 접근법인 입자법은 격자 생성이 불필요하며, 이류항(Convection)의 계산을 입자의 이동으로 직접 계산함으로써 격자법에서 심각하게 유발되는 수치 확산이 없을 뿐 아니라, 질량이 완벽히 보존된다는 장점을 갖고 있다. 또한 격자 생성에 소요되는 시간과 노력을 최소화 할 수 있다. 이러한 입자법의 예로는 크게 SPH(Smooth particle hydrodynamics)법(Monaghan, 1988)과 MPS(Moving particle simulation)법(Koshizuka and Oka, 1996)이 있으며, 격심한 자유표면 운동(Marrone et al., 2012), 혼상유동(Khayyer and Gotoh, 2013; Jeong et al., 2013), 유체-고체 연성 해석(Yang et al., 2012) 등의 연구에 활발히 응용되고 있다.

하지만, 상술한 입자법의 경우 SPH법이나 MPS법 모두 압력의 시간 및 공간적 진동에 의해 유동장이 불안해 지는 문제점들이 지적되고 있다. 이러한 극심한 압력 진동은 비물리적이며 유체의 거동을 부자연스럽게 표현할 가능성이 있다. 따라서 유동장의 타당한 해석을 위해서는 압력 진동의 안정화가 중요한 과제이다. 압력 진동을 유발시키는 가장 큰 원인 중 하나는 유체 입자가 외력에 의해 운동할 경우 매 시각 한 입자의 유효반경 안에 포함되는 주변입자들의 공간적인 분포가 불규칙해짐에 따라 지배방정식의 편미분 연산자에 대한 기존의 근사 모델들에 대한 정확성이 보장되지 않는다는 점에 있다. 이를 해결하기 위하여, 특히 비압축성 유동해석을 하기 위한 MPS법의 경우, 많은 연구자들에 의한 다양한 연구결과가 제시되었다(Tanaka and Masunaga, 2008; Lee et al, 2011; Khayyer and Gotoh, 2011; Khayyer and Gotoh, 2012). 이 가운데 Koh et al.(2012)은 지배방정식의 각 편미분항을 Taylor 급수 전개를 통한 이동최소자승(MLS, Moving least square) 근사법으로 이산화하여 입자의 불규칙한 분포로 인해 발생하는 공간적인 수치오차를 줄 일 수 있는 CPM(Consistent particle method)법을 제안하였다.

본 논문에서는 Koh et al.(2012)가 제안한 CPM법의 MLS 근사에 비해 보다 가까운 위치에 있는 입자정보에 가중치를 부여하는 가중이동최소자승(WMLS, Weighted moving least-squares) 기반의 입자법 시뮬레이션 기술을 개발한다. 일반적으로 WMLS는 MLS에 비해 유효반경 내의 보다 가까운 위치에 있는 물리량에 큰 가중치를 주게 되므로 보다 효과적으로 예측이 가능하다는 장점이 있다. 개발된 시뮬레이션 기술의 타당성을 살펴보기 위하여 WMLS법을 이용해 구한 편미분항의 연산자를 2차원의 해석해를 갖는 간단한 문제를 통해 검증하였으며, 정수압 문제 및 댐 붕괴 문제에 적용하여 공학적 적용성을 검토하였다.

2. 수치 시뮬레이션 기술

2.1 지배 방정식

비압축성 및 점성유동에 관한 지배 방정식은 다음과 같은 연속방정식과 Navier-Stokes 방정식이다.

여기서 ρ는 밀도, t는 시간, u는 속도벡터, ∇는 구배 연산자, P 는 압력, v는 동점성계수, f는 외력을 각각 나타낸다.

식 (2)의 좌변은 이류항을 포함한 라그란지의 전미분(Total derivative) 형태이고, 이는 입자의 직접적인 이동에 의해 계산된다. 한편, 우변은 압력 구배항, 점성항 그리고 외력항으로 구성되어 있다. 지배방정식의 편미분 연산자를 포함하는 모든 항은 후술하는 WMLS 근사를 이용하여 이산화한다.

2.2 WMLS에 기반한 편미분 연산자 근사

중심입자 i(xi,yi) 주변의 물리량 분포를 Taylor 급수 전개를 하면 식 (3)과 같다.

단, h = xj–xi , k = yj–yi , 그리고 우변에서 물리량의 공간미분은 다음과 같이 정의된다.

식 (3)에서 주변 입자들을 고려하면 다음의 연립방정식이 유도된다.

단, N 은 영향 반경 내의 주변입자의 총 개수이다.

최소자승(LS, Least-squares) 근사를 비롯한 각종 근사 추정법에 있어서, 보간된 정보의 취득방법에 따라 계산이 불안정해져 수치적 진동을 발생시키는 경우도 있다. 따라서 추정오차를 줄이기 위해 가급적 중심 입자와 가까운 주변 입자들을 이용하는 것과 동시에, 유동의 정보를 모든 방향에서 얻는 것이 바람직하다. 본 연구에서는, 거리에 따라 가중치를 부여하는 WMLS 근사법을 채택한다(Jeong et al., 2009). 가중함수 wj를 적용한 잔차의 식 (8)과 최소화 조건인 식 (9)로부터 식 (10)과 같은 연립방정식이 유도된다. 단, A 는 계수행렬로서 식 (11)과 같이 주어진다.

최종적으로, 편미분 연산자를 구하기 위해서는 식 (10)의 계수행렬 A의 역행렬 계산이 필요하며, 본 연구에서는 QR-해법을 이용하였다. 또한, 행렬 내의 제 1, 2행은 구배 연산자, 제 3, 5행은 확산 연산자에 관여한다.

본 연구에서 사용된 가중함수는 PNU-MPS(Pusan-National-University-modified Moving Particle Simulation)법(Lee et al., 2011)과 동일하다.

2.3 비압축성 알고리즘

만약 입자의 질량이 모두 동일하다고 가정하면, 유체의 밀도는 입자수밀도에 비례하므로 연속방정식 식 (1)은 입자수밀도가 일정하다는 조건과 동일하다.

본 방법에서는 비압축성 유동의 계산 알고리즘으로서 PNU-MPS법(Lee et al., 2011)과 유사한 알고리즘을 사용하며, Fig. 1에는 계산 과정을 나타낸다. 매 시간 스텝의 초기에는 편미분 연산자의 계수를 계산하며, 그 후 연산은 양과 음의 2단계로 나눈다. 시각 n에 있어서 입자의 위치, 속도, 압력을 각각 , , 이라고 하면, 제 1단계에서는 점성항과 외력항의 양적 계산을 통하여 입자의 중간 속도 를 계산한다.

Fig. 1

Procedure of numerical algorithm

다음으로, (Poisson) 방정식을 음적으로 계산한다.

단, 위 식의 우변의 상수 𝛾 는 비압축성 조건을 만족하는 범위내에서 1.0 이하의 값이다.

식 (13)의 좌변은 확산 모델에 의해 연립 1차방정식으로 이산화 할 수 있으며, Jeong et al.(2009)에서 제안한 SOR(Successive over relaxation)에 의한 반복해법에 의해 구하였다. 단, SOR 반복해법에 사용한 완화계수(Relaxation factor)는 Jeong et al. (2009)에서 사용한 1.5를 사용하였다.

압력이 구해진 뒤 식 (14)로부터 속도의 수정치 를 계산한다. 단, 식 (14)의 우변의 압력 구배 계산은 구배모델을 사용한다.

최종적으로, 시각 n+1 에서 입자의 속도와 위치는 각각 다음 두 식에 의해 새롭게 얻어진다.

Fig. 2

Arc method

2.4 자유표면 경계조건

자유표면 경계조건으로 운동학적 경계조건과 동역학적 경계조건이 있다. 본 연구에서 운동학적 경계조건은 자유표면 입자를 직접 이동하여 만족시키고, 동역학적 경계조건은 자유표면 입자의 압력을 대기압과 동일한 0의 값으로 고정함으로써 만족시킨다. 두 경계조건을 만족시키기 위해서는 자유표면 입자의 탐색이 선행되어야 한다. 본 연구에서는 원의 기학학적인 특징을 이용하여 Dilts(2000)가 제안한 “Arc”의 원리를 이용한 자유표면 입자탐색법(Park et al., 2014)을 적용하였다. 우선, Fig. 5와 같이 초기 입자간 거리 보다 큰 직경 R 을 각 입자마다 설정한다. 이 때 입자 A 의 경우 주변 입자들에 의해 원의 모든 호가 닫히게 되어 유체 내부 입자로 판단되는 반면, 입자 B 의 경우 원호의 일부가 닫히지 않아 자유표면 입자로 판단된다. 본 논문에서 직경 R 은 Park et al.(2014)과 동일하게 2.1l0를 사용하였다.

Fig. 5

Comparison of Laplacian model

3. 수치 시뮬레이션

3.1 구배 및 확산 모델 검증

본 연구에서 개발한 구배 모델과 확산 모델의 타당성을 확인하기 위해 2차원의 간단한 문제에 적용시켜 보았으며, 본 계산 알고리즘과 유사한 PNU-MPS법(Lee et al., 2011)의 구배 및 확산 모델과 비교하였다. 단, PNU-MPS법의 구배 및 확산 모델은 다음과 같다.

여기서 d는 계산 공간의 차원, λ 는 해석해와 동일한 분산을 유지하기 위해 사용하는 계수이다.

Fig. 3과 같은 계산 영역 [0.8, 4.1]⨉[0.8, 4.1]에서 입자배치를 정규 입자배치의 상태 기준으로 50%의 무질서도(Randomness) 를 주었다. 초기 물리량은 식 (20)과 같이 설정하였으며, 구배모델과 확산 모델의 해석해는 각각 식 (21), (22)와 같다. 초기 입자간 거리 l0 는 0.1m이며, 임계거리 re는 3·1l0, 경계조건은 경계면에서 물리량의 기울기가 영이라는 Neumann 조건을 적용하였다.

Fig. 3

Example of irregular nodes with random ratio of 50%

Figs. 4~5는 각각 해석해, PNU-MPS법 그리고 본 연구의 구배 및 확산 모델을 적용한 결과이며, Figs. 6~7은 수치해석결과와 해석해와의 물리량의 차이를 나타낸다. 이 때 구배모델의 결과는 x, y방향의 총합을 나타낸다. 전체적으로 PNU-MPS법의 경우 공간적인 진동을 보이며 해석해와의 차이를 보이는 반면, 본 연구의 구배 모델과 확산 모델의 경우 해석해와 잘 일치하는 것을 확인 할 수 있다.

Fig. 4

Comparison of gradient model

Fig. 6

Difference of gradient values between analytic and numerical results

Fig. 7

Difference of Laplacian values between analytic and numerical results

다음으로, 보다 정량적인 해석을 위하여 다음 식을 이용하여 해석해에 대한 상대오차율을 분석하였다.

PNU-MPS법의 경우 구배 모델과 확산 모델의 상대오차율은 각각 68.5%와 60.4%인. 반면, 본 연구의 경우 상대오차율이 각 각 6.5%와 3.2% 정도로 해석해와 잘 일치하는 경향을 보인다. 따라서, 본 연구의 구배 및 확산 모델은 입자의 불규칙한 배열에서도 정확도가 높게 계산된다는 것을 확인 할 수 있다.

3.2 정수압 시뮬레이션

본 연구에서 개발된 방법을 이용하여 Fig. 8에 보이는 바와 같이, 수조의 폭과 수심이 각각 0.4m 인 2차원 사각 수조내의 정수압 문제에 적용하였다. 사용한 총 입자 개수는 약 1700개이며 이 중 물 입자는 1600개이다. 중력가속도와 밀도는 각각 9.81m/s2과 1000kg/m3으로 설정하였다. 총 계산시간은 5.0초이다. 비교를 위하여 수조 바닥의 중심 P1 지점에서 압력을 측정하였다.

Fig. 8

Schematic description of initial setup for hydo-static pressure test

Fig. 9는 P1 지점에서 계측한 압력의 시간변화를 나타낸다. PNU-MPS법의 경우 고주파수의 비물리적인 압력 진동을 포함하고 있으며, 특히 1초 이후에 해석해와 다소 차이를 보인다. 반면 본 시뮬레이션의 경우 유동장의 압력이 안정이 되어 진동이 거의 없으며 정확도도 크게 개선되어 해석해와 잘 일치하는 것을 알 수 있다.

Fig. 9

Time history of pressure probed at P1

3.3 댐붕괴 시뮬레이션

댐붕괴 문제는 자유표면 거동에 비선형성 특성을 포함하므로 비선형성 자유표면 유동을 포함하는 다양한 수치 시뮬레이션의 검증에 이용되어 왔다. 따라서, 본 절에서는 Martin and Moyce(1985)가 수행한 2차원 댐붕괴의 조건과 동일한 조건하에 서 시뮬레이션을 수행하였다. Fig. 10과 같이 초기 형상을 갖는 2차원의 수조 내부의 왼쪽에는 폭 0.4m, 높이 0.8m의 물기둥이 위치한다. 점성의 영향은 고려하였고 표면장력의 영향은 무시하였다. 시뮬레이션의 결과 비교를 위하여 수조의 바닥으로부터 0.02m 떨어진 우측 벽면의 P2 지점에서 압력을 계측하였다.

Fig. 10

Schematic description of initial setup for dam-broken simulation

Fig. 11은 P2 지점에서 시간에 따른 압력변화를 나타낸다. 시뮬레이션 결과, 시간에 따른 압력의 경향은 PNU-MPS법의 결과와 유사하게 나타나지만, 본 시뮬레이션의 경우, 유체가 벽면에 충돌하는 압력의 첫 번째 극치(Peak) 값 부근의 압력 진동이 PNU-MPS법에 비해 안정됨을 알 수 있다. Fig. 12는 붕괴된 수괴가 오른쪽 벽면으로부터 반사되어 다시 붕괴되는 순간의 모습을 비교하여 나타낸다. 전반적으로 붕괴파의 첨단(Front) 형상을 포함한 자유표면 형상은 유사해 보이지만 본 시뮬레이션의 경우 PNU-MPS법에 비해 첨단 입자들이 보다 집중적인 것을 볼 수 있다. 또한 오른쪽 벽면 위로 치솟는 산발적인 비말 (Splash)의 표현이 완화되어 나타나 있다. 하지만, 몇몇 입자들의 거센 반발력에 의해 공중에 흩어지는 것을 볼 수 있는데, 이는 비물리적인 것으로 불 수 있으며 자유표면상의 입자들의 충돌모델의 개선이 필요한 것으로 판단된다.

Fig. 11

Time history of pressure probe at P2

Fig. 12

Comparison of wave configuration at t=1.26s

Fig. 13은 붕괴되는 물기둥 첨단의 시간에 따른 위치변화를 실험(Martin and Moyce, 1985), 타 시뮬레이션인 SOLA-VOF법 (Hirt and Nicols, 1981)과 PNU-MPS법(Lee et al., 2011)와 비교하여 나타낸다. 단, 실험의 경우 1.125와 2.5는 초기 물기둥 크기의 종횡비를 나타낸다. 본 시뮬레이션의 경우 붕괴파의 첨단 (Front) 부분이 수조 바닥에서 수평방향으로의 진행속도가 PNU-MPS법에 비해 다소 지연됨이 관찰되었으며, Martin and Moyce(1985)의 실험에 보다 근접한 결과를 보여주고 있다.

Fig. 13

Comparison for position change of leading edge of collapsed water column

4. 결 론

본 논문에서는 가중이동최소자승(WMLS, Weighted moving Least-squares)법에 기반한 입자법 시뮬레이션 기술을 개발하였다. 이는, Koh et al.(2012)가 제안한 CPM법의 이동최소자승 (MLS) 근사에 비해 유효반경 내의 보다 가까운 위치에 있는 물리량에 큰 가중치를 주게 되므로 보다 효과적으로 예측이 가능하다는 장점이 있다. 따라서 지배방정식의 각 편미분항을 WMLS법을 이용하여 이산화 하였다. 본 방법의 검증을 위해 WMLS법을 통해 구한 편미분항의 연산자인 구배 및 확산 모델을 2차원의 해석해를 갖는 단순한 문제에 적용하여 그 타당성을 확인하였다. 또한, 정수압 문제와 댐붕괴 시뮬레이션을 수행하여, 해석해, 실험 및 타 시뮬레이션 결과와 비교하였다. 결과적으로 본 연구에서 개발한 입자법의 경우 유동장의 압력장이 크게 안정적이며, 정확도가 개선되었음을 알 수 있었다. 향후 본 연구에서 개발한 수치 스킴의 유용성을 확인하기 위한 기초적인 수치실험과, 해양공학 분야의 보다 다양한 응용문제에 적용하여 그 타당성을 검증할 필요가 있을 것이다.

Acknowledgements

본 연구는 2011년도 정부(교육과학기술부)의 재원으로 한국연구재단의 기초연구사업 지원을 받아 수행된 것(2011-0011220)을 밝히며, 연구비 지원에 감사드립니다.

References

Dilts G.A.. Moving Least-squares Particle Hydrodynamics II: Conservation and Boundaries. International Journal for Numerical Methods in Engineering 48(10)2000;:1503–1524.
Hirt C.W., Nichols B.D.. Volume of Fluid(VOF) Method for the Dynamics of Free Boundaries. Journal of Computational Physics 391981;:201–225. 10.1016/0021-9991(81)90145-5.
Jeong S.-M., Nam J.-W., Hwang S.-C., Park J.-C., Kim M.-H.. Numerical Prediction of Oil Amount Leaked from a Damaged Tank using Two-dimensional Moving Particle Simulation Method. Ocean Engineering 692013;:70–78. 10.1016/j.oceaneng.2013.05.009.
Jeong S.M., Park J.C., Heo J.K.. Numerical Study on Two-Dimensional Incompressible Viscous Flow based on Gridless Method. Journal of Computational Fluids Engineering 14(4)2009;:93–100.
Koh C.G., Gao M., Luo C.. A New Particle Method for Simulation of Incompressible Free Surface Flow Problems. International Journal for Numerical Methods in Engineering 89(12)2012;:1582–1604. 10.1002/nme.3303.
Koshizuka S., Oka Y.. Moving-particle Semi-implicit Method for Fragmentation of Incompressible Fluid. Nuclear Science and Engineering 1231996;:421–434.
Khayyer A., Gotoh H.. Enhancement of Stability and Accuracy of the Moving Particle Semi-implicit Method. Journal of Computational Physics 2302011;:3093–3118. 10.1016/j.jcp.2011.01.009.
Khayyer A., Gotoh H.. A 3D Higher Order Laplacian Model for Enhancement and Stabilization of Pressure Calculation in 3D MPS-based Simulations. Applied Ocean Research 372012;:120–126. 10.1016/j.apor.2012.05.003.
Khyyer A., Gotoh H.. Enhancement of Performance and Stability of MPS Mesh-free Particle Method for Multiphase Flows Characterized by High Density Ratios. Journal of Computational Physics 2422013;:211–233. 10.1016/j.jcp.2013.02.002.
Lee B.-H., Park J.-C., Kim M.-H., Hwang S.-C.. Step-by-step Improvement of MPS Method in Simulating Violent Free-surface Motions and Impact-loads. Computer Methods in Applied Mechanics and Engineering 2002011;:1113–1125. 10.1016/j.cma.2010.12.001.
Marrone S., Bouscasse B., Colagrossi A., Antuono M.. Study of Ship Wave Breaking Patterns using 3D Parallel SPH Simulations. Computers & Fluids 692012;:54–66. 10.1016/j.compfluid.2012.08.008.
Martin J.C., Moyce W.J.. An Experimental Study of the Collapse of Liquid Columns on a Rigid Horizontal Plane. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 2441985;:312–324.
Miyata H., Park J.C., Rahman. M.. Ch.5 Wave Breaking Simulation. Potential Flow of Fluids Computational Mechanics Publications. UK: 1995. p. 149–176.
Monaghan J.J.. An Introduction to SPH. Computer Physics Communications 481988;:89–96. 10.1016/0010-4655(88)90026-4.
Park J.-I., Park J.-C., Hwang S.-C., Heo J.-K.. Two-Dimensional Particle Simulation for Behaviors of Floating Body near Quaywall during Tsunami. Journal of Ocean Engineering and Technology 28(1)2014;:12–19. 10.5574/KSOE.2014.28.1.012.
Sussman M., Smereka P., Osher S.. A Level Set Approach for Computing Solutions to Incompressible Two-phase Flow. Journal Computational Physics 1141994;:272–280.
Takewaki H., Yebe T.. The Cubic-Interpolated Pseudo Particle(CIP) Method: Application to Nonlinear and Multidimensional Hyperbolic Equations. Journal of Computational Physics 70(2)1987;:355–372. 10.1016/0021-9991(87)90187-2.
Tanaka M., Masunaga T.. Stabilization Smoothing of Pressure on MPS Method by Quasi-compressibility. Transactions of the Japan Society for Computational Engineering and Science 20082008;:20080025.
Yang Q., Jones V., McCue L.. Free-surface Flow Interactions with Deformable Structures using an SPH. FEM model. Ocean Engineering 552012;:136–147. 10.1016/j.oceaneng.2012.06.031.

Article information Continued

Fig. 1

Procedure of numerical algorithm

Fig. 2

Arc method

Fig. 5

Comparison of Laplacian model

Fig. 3

Example of irregular nodes with random ratio of 50%

Fig. 4

Comparison of gradient model

Fig. 6

Difference of gradient values between analytic and numerical results

Fig. 7

Difference of Laplacian values between analytic and numerical results

Fig. 8

Schematic description of initial setup for hydo-static pressure test

Fig. 9

Time history of pressure probed at P1

Fig. 10

Schematic description of initial setup for dam-broken simulation

Fig. 11

Time history of pressure probe at P2

Fig. 12

Comparison of wave configuration at t=1.26s

Fig. 13

Comparison for position change of leading edge of collapsed water column