효율적인 해양구조물 유동 해석을 위한 직교좌표계 기반의 코드 개발 - AMR, VOF, IBM, VIV, LES의 통합

doi:10.3744/SNAK.2014.51.5.409

OA학술지
Journal of the Society of Naval Architects of Korea

효율적인 해양구조물 유동 해석을 위한 직교좌표계 기반의 코드 개발 - AMR, VOF, IBM, VIV, LES의 통합 Development of a Cartesian-based Code for Effective Simulation of Flow Around a Marine Structure - Integration of AMR, VOF, IBM, VIV, LES

DOI : 10.3744/SNAK.2014.51.5.409
Publish: Journal of the Society of Naval Architects of Korea Volume 51, Issue5, p409~418, 20 Oct 2014

ABSTRACT

Simulation of flow past a complex marine structure requires a fine resolution in the vicinity of the structure, whereas a coarse resolution is enough far away from it. Therefore, a lot of grid cells may be wasted, when a simple Cartesian grid system is used for an Immersed Boundary Method (IBM). To alleviate this problems while maintaining the Cartesian frame work, we adopted an Adaptive Mesh Refinement (AMR) scheme where the grid system dynamically and locally refines as needed. In this study, We implemented a moving IBM and an AMR technique in our basic 3D incompressible Navier-Stokes solver. A Volume Of Fluid (VOF) method was used to effectively treat the free surface, and a recently developed Lagrangian Dynamic Subgrid-scale Model (LDSM) was incorporated in the code for accurate turbulence modeling. To capture vortex induced vibration accurately, the equation for the structure movement and the governing equations for fluid flow were solved at the same time implicitly. Also, We have developed an interface by using AutoLISP, which can properly distribute marker particles for IBM, compute the geometrical information of the object, and transfer it to the solver for the main simulation. To verify our numerical methodology, our results were compared with other authors' numerical and experimental results for the benchmark problems, revealing excellent agreement. Using the verified code, we investigated the following cases. (1) simulating flow around a floating sphere. (2) simulating flow past a marine structure.

KEYWORD

가상경계법 , 유체체적법 , 적응격자기법 , 대와류모사 , 와류유발진동 , 해양구조물

본문

Collapse all

1. 서 론

해양구조물 주위 유동해석에 있어 기존 해석 코드들은 크게 선형 포텐셜(potential) 기반의 해석 코드와 상용코드(commercial code)로 구분할 수 있다. 하지만 이런 코드들은 다상 유체와 구조물 사이의 상호작용 해석에 있어 한계를 지닌다. 선형 포텐셜 기반의 해석 프로그램들은 수치 감쇠가 적어 생성된 파를 멀리 전파할 수 있다는 장점이 있지만 자유표면이 쇄파하거나 분리, 결합하는 비선형 현상을 표현할 수 없고, 구조물 근처에서 점성효과를 정확히 반영할 수 없다. 또한 슬래밍(slamming) 현상 및 높은 파고의 파도를 구현할 수 없으며 난류 해석이 불가능하다. 상용코드를 이용한 해석은 일반적으로 격자생성이 매우 복잡할 뿐만 아니라, 해석 결과의 정밀도와 전산 효율이 저하되는 문제점이 있다. 이를 완화시키기 위하여 비정렬 격자계(unstructured grid system)가 사용되고 있으나, 격자점 배치의 효율성은 증대된 반면 구현의 복잡성, 다른 코드와의 호환성 저하, 계산의 수렴성 저하 등의 문제점은 여전히 남아 있다. 또한 해양구조물 관련 점성 유동인 경우 대부분이 높은 레이놀즈 수(reynolds number)에서 해석해야 하는 어려움이 있다. 높은 레이놀즈 수에 의한 난류를 정확히 예측하기 위해서는 많은 격자점과 정밀한 LES 난류모델을 필요로 한다.

본 연구에서는 복잡한 형상의 해양구조물을 직교좌표계에서 구현하기 위해 IBM(Immersed Boundary Method) 기법을 사용하였다. 기존 IBM에서는 복잡한 형상의 물체를 구현하기 위해서 물체 표면을 나타내는 방정식을 이용하였지만, 이런 방법의 경우 간단한 형상에만 적용할 수 있으며, 형상이 복잡해지면 정확한 물체의 표면을 구현하지 못한다는 어려움이 있다. 이러한 문제를 해결하기 위해서 본 연구에서는 표면점(marker particle)을 이용하는 방법을 채택하였다 (Choi, et al., 2011). 이 방법은 3차원 모델링 상용프로그램인 AutoCAD로 모델링된 형상의 표면점들을 추출하여 Solver에 전달하고, 전달된 표면점들의 좌표정보를 이용하여 복잡한 형상을 구현할 수 있다. 또한 형상이 변경되어도 Solver의 수정 없이 AutoCAD 모델링만 변경하면 해석이 가능하다.

표면점 기술이 적용된 IBM과 함께 효율적인 격자 구조 기법인 AMR(Adaptive Mesh Refinement) 기법을 사용하였다. 복잡한 형상의 구조물 주위에는 많은 격자점을 필요로 하며, 구조물의 운동을 포함하여 해석할 경우 구조물의 움직임에 따라서 조밀한 격자를 사용해야 하기 때문에 전체적인 컴퓨터 메모리양과 수치계산 시간의 급격한 상승을 초래한다. 또한 다상유동의 대표적 특징인 자유표면을 정확히 구현하기 위해서는 자유표면 근처에서 조밀한 격자를 사용해야 한다. 이를 극복하기 위해 사용된 AMR기법은 격자구조에 단계(level)를 두어 물리량(속도, 압력)들의 변화가 심하지 않은 구역에서는 성긴 격자를 사용함으로써 컴퓨터 메모리 및 계산시간을 절감할 수 있으며, 물리량들이 급변하는 곳에서는 조밀한 격자를 사용함으로써 더욱 정확한 계산을 수행할 수 있다.

AMR 격자구조에서 출발하여 표면점 기술이 적용된 IBM 및 높은 수의 난류를 정확히 해석하기 위한 LES 해석 기술을 통합하였다. 또한 VOF(Volume of Fluid) 기법을 적용하여 자유표면을 포함하는 다상 유체와 구조물과의 상호작용을 정확히 해석할 수는 코드를 완성하였다. 완성된 코드를 이용하여 자유표면에 떠 있는 구의 동적 거동 및 수면의 변화에 대한 연구와, 부분적으로 잠수된 해양구조물 주위의 수면의 변화 및 구조물이 받는 힘에 대한 연구를 수행하였다.

2. 수식해석 기법

   2.1 지배방정식

구조물의 움직임이 포함된 다상유동 해석을 위한 가상경계법이 적용된 LES의 지배 방정식은 연속 방정식, 운동량 방정식, 이송방정식, 구조물의 동역학 방정식으로 구성되며, 다음과 같다.

여기서 는 box filter를 사용하여 여과된 속도 성분이고, 는 압력, 𝞽_ij는 레이놀즈 응력이다. 운동량 방정식에서 는 가상경계법에 의해 추가된 운동량 부가, F_b는 체적 힘을 의미한다. 체적힘은 중력에 의한 힘이며, 중력가속도는 9.81m/s이다. 구조물의 운동방정식은 1차 상미분방정식 형태로 표현할 수 있다(식 (4)). 운동 방정식에서 q₁ = [ q₁ q₂ ··· q_n ]^T은 동역학 시스템에서의 n개의 일반화된 좌표계 집합을 의미하며, 는 일반화된 속도 집합을 의미한다. I 는 n × n 단위행렬(identity matrix), 0은 n × n 영행렬이다. M (q₁) 은 질량행렬, F_f(q₁, q₂,t) 는 운동체에 가해지는 힘을 의미하며, 본 연구에서는 중력과 유체력만 고려하였다.

각 방정식은 유한체적법(finite volume method)으로 차분되며, 비균일 엇갈린 격자계(nonuniform staggered grid)를 사용하였다. 연속방정식과 운동량 방정식을 분리하기 위하여 Fractional Step Method (Kim & Moin, 1985)을 사용하였다.

   2.2 Adaptive Mesh Refinement (AMR)

본 연구에 사용된 AMR 기법 (Vanella, et al., 2010)에서 계산영역은 블럭으로 구성되어 있으며 각 블럭은 nx × n𝓎 × nz개의 균일한(uniform) 격자 구조를 가진다(Fig. 1). 블럭의 단계(level) Ω^l은 격자 구조의 크기를 결정한다. 즉 가장 성긴 격자 구조를 가지는 블럭의 단계는 Ω¹이며, 가장 조밀한 격자 구조를 가지는 블럭의 단계는 이다. 임의의 Ω^l 블럭이 Ω^l+1 블럭으로 분리될 때 Ω^l 블럭을 parent 블럭, Ω^l+1 블럭을 children 블럭이라 한다. Children 블럭들의 합과 parent 블럭은 정확히 같은 영역을 차지한다. 각각의 블럭에서, 각 방향으로의 셀(cell)의 크기는 이다. 이때 아래첨자 i는 방향을 나타내며, 위첨자 l은 블럭의 단계를 나타낸다. AMR 기법에서 가장 중요한 요소는 단계가 다른 블럭들이 인접해 있을 때 속도와 압력 같은 변수들이 경계에서 잘 전달되는지의 여부이다. 본 연구에서는 Restriction과 Prolongation 알고리즘 (Vanella, et al., 2010)을 사용하여 단계가 다른 블럭에서 경계값을 정확히 계산하여 물리량들이 정확히 전달될 수 있게 하였다.

[Fig. 1] AMR block system

   2.3 자유표면

본 연구에서는 자유표면을 해석할 수 있는 방법 중 하나인 VOF기법을 사용하여 자유표면을 구현하였다. VOF 방법은 각 cell마다 정의되는 용적함수 f(volume fraction)를 시간에 따라 추적하여 경계면을 재구성하는 방법으로, 계산과정에서 질량이 정확히 보존된다는 장점이 있어 가장 보편적으로 사용되어지고 있다. 설명의 편의를 위해서 두 물질을 액체, 기체 즉 이상유동(Two-Phase flow)으로 가정하면, f는 각 cell에서 차지하는 액체의 부피비로 다음과 같이 정의한다.

즉 f값이 1이면 액체, 0이면 기체, 0과 1 사이면 두 물체의 경계를 나타낸다. 각 cell에서 물질들의 부피비가 정해지면, 밀도와 점성계수는 다음과 같이 정의된다.

이때 ρ₁, μ₁는 액체의 밀도와 점성계수, ρ₂, μ₂는 기체의 밀도와 점성계수이다.

VOF 방법에서는 실제로 경계면을 추적하는 것이 아니라 용적함수 f값을 추적하게 된다. 특정 시간에 공간에 대하여 주어진 스칼라 함수인 f값을 이용하여 경계면을 형성하는 과정이 필요한데 본 연구에서는 ELVIRA(Effective Least squares Volume-of-fluid Interface Reconstruction Algorithm; Pilliod, 1992)을 사용하였다. 또한 이송방정식의 차분은 Operator split advection algorithm (Pilliod & Puckett, 2004)을 사용하였다.

   2.4 와류유발진동 (VIV)

구조물의 운동 방정식과 유체 지배 방정식이 서로 연동되어 있기 때문에 VIV 해석은 쉽지 않다. 가장 쉽게 연동된 방정식들의 해를 구하는 방법은 물체의 운동 방정식을 계산하기 위해서 필요한 힘을 시간에 대해 양해적(explicit)으로 분리하는 것이다 (Hall, 1994; Robertson, et al., 2003; Jadic, et al., 1998). 즉, 각 시간 스텝(step)에서 유체지배 방정식으로 부터 계산된 힘을 운동방정식에 대입하여 물체의 위치와 속도를 결정해주는 방법이다. 그러나 이러한 방법은 물체와 유체 사이에 발생하는 상호작용을 완벽히 구현해내지 못한다. 정확한 VIV 해석을 위해서는 물체의 운동 방정식과 유체지배 방정식을 시간에 대해 음해적(implicit)으로 동시에 계산해야 한다. 이 과정에서 각 시간 전진마다 반복 해석하는 기술이 요구되며, 이로 인해 많은 계산 시간을 필요로 하나, 상호작용을 보다 완벽히 구현해낼 수 있다는 장점을 가지고 있다. 본 연구에서는 Hamming`s 4th-order, predictor-corrector method (Yang & Balaras, 2006)을 이용하여 음해적으로 계산을 수행하였다.

   2.5 난류모델

해양구조물과 같이 구조물의 크기가 큰 경우에는 대부분이 높은 레이놀즈수의 난류 특성을 보인다. 따라서 이런 난류유동을 시간에 따라 정확히 해석할 수 있는 LES 기법이 필요하다. 본 연구에서는 Lagrangian Dynamic Subgrid-scale Model (LDSM; Meneveau, et al., 1996)모델을 사용하였다.

3. 수치해석 기법의 검증

   3.1 구의 자유낙하

수치해석 코드의 검증을 위하여 자유 낙하하는 구를 대상으로 종단 속도와 종단 항력을 이론값과 비교해 보았다. 계산 영역은 Fig. 2와 같다. 낮은 속도(Re≪1)에서 구에 작용하는 항력(= 3πμDV )은 Stokes Law를 통해 계산할 수 있으며, 이 힘은 구의 무게 W(=ρ_s × (πD³/6) × g)와 부력 F_b(=ρ_f × (πD³/6) × g)과 힘의 균형을 이룬다(식 (8)).

[Fig. 2] Physical configuration and computational domain of the free-falling sphere problem

이때 D는 구의 지름, V는 자유 낙하 속도, ɡ는 중력가속도의 크기, ρ_s와 ρ_f는 각각 구와 유체의 밀도이다. 점성계수(μ) 100kg/(ms)인 가상의 유체 속에서 밀도(ρ_s) 171.0 kg/m³인 구를 자유 낙하하면서 구의 속도와 항력을 계산해 보았다.

계산 영역은 구의 지름(1m)을 기준으로 하여 x, 𝓎, z방향으로 각각 80m, 20m, 20m이며 초기 구 중심의 위치는 각 방향으로 10m이다. 중력은 +x 방향으로 항력과 부력은 −x 방향으로 설정하였다. 모든 면에서는 Dirichlet조건(u=𝜐=w=0)을 사용하였다. 구 주위에는 조밀한 블럭을 배치하였고, 구와 멀리 떨어진 구역에는 성긴 블럭을 배치하여 총 480개의 블럭을 사용하였다. 각 블럭은 각 방향으로 8개의 격자들로 구성된다. 본 계산에 사용된 dt는 0.0001초이다.

식 (8)으로부터 계산된 종단속도는 0.9262m/s이다. 종단 속도를 기준으로 계산한 Re(=ρ_fV_terminalD/μ)는 9.262×10E−3으로 Stokes law의 가정인 Re≪1에 잘 부합된다.

Fig. 3은 자유 낙하하는 구의 속도와 항력을 나타낸 그림이다. 종단 속도의 99%에 도달하는 이론적인 시간은 0.44초 이며, 그때의 이론 속도는 0.917m/s이다. 본 연구에서는 0.44초에서의 구의 속도가 0.905m/s로 약 1.3%의 오차로 이론값과 잘 일치하는 것을 확인할 수 있다. 또한 항력의 크기도 이론값과 잘 일치하는 것을 확인할 수 있다(Fig. 3(b)).

[Fig. 3] Results of the free-falling sphere problem: (a) time history of the velocity, (b) time history of the drag coefficient

   3.2 Broken water-dam

두 번째 검증으로 가상의 댐이 갑자기 사라졌을 경우, 가로 1m, 세로 2m인 물기둥의 거동을 예측하는 문제이다. 전체 계산 영역은 4m × 4m이며, 물성치는 공기의 경우 밀도 1.225 kg/m³, 점성계수는 1.7763×10⁻⁵kg/ms이고, 물의 경우 밀도는 999.2kg/m³, 점성계수는 1.1377×10⁻³kg/ms을 사용하였다(Fig. 4). 계산에 사용된 dt는 10⁻⁶이다.

[Fig. 4] Physical configuration of a broken water-dam problem

Fig. 5는 물기둥이 무너져 내릴 때 물기둥의 최대 변위 Z를 초기값으로 나눈 것을 시간에 대하여 나타낸 것이다. 다른 연구자들의 결과 (Park, et al.., 2003; Martin & Moyce, 1952; Koshizuka & Oka, 1996; Hirt & Nichols, 1981)와 잘 일치하는 것을 알 수 있다.

[Fig. 5] Comparison of length of the surge front(z/？) against time

4. 다상유체와 구와의 VIV 해석

자유표면에서 자유롭게 운동하는 구의 VIV 해석을 수행하였다. 구의 움직임에 따라 자유표면의 모양 및 유체력이 변하게 되고 이는 다시 구의 움직임에 영향을 미치게 된다. 다상유체와 구의 VIV 거동 해석은 단상에서의 VIV 문제보다 난이도가 훨씬 높다. 단상에 비해 추가적으로 자유표면을 정확히 구현할 수 있어야 한다. 또한 구의 거동을 정확히 해석하기 위해서는 구가 받는 힘을 정확하게 계산하여야 한다.

   4.1 계산영역 및 경계조건

전체 계산영역 및 경계조건은 Fig. 6에 나타내었다. 전체 계산영역의 크기는 각 방향으로 3m×2m×2m이며 구의 지름 D는 0.4m이다. 바닥 면(−z방향 면)에는 no-slip조건을 사용하였고, 나머지 면에는 slip 경계조건을 사용하였다. 자유표면의 초기 위치는 바닥면에서 1.2m 떨어져 있으며, 이를 z축 원점으로 설정하였다. 초기 구는 정지 상태이며 자유표면에 정확히 반이 잠겨있다. 계산에 사용된 물의 물성치는 밀도 ρ_water =1000kg/m³, 점성계수 μ_water =1.141× 10⁻³kg/ms이며, 공기의 경우 밀도 ρ_air =1.2kg/m³, 점성계수 μ_air =1.791× 10⁻⁵kg/ms을 사용하였다. 구의 밀도는 물의 밀도보다 가벼운 ρ_sphere =900kg/m³이다. 중력가속도는 −z방향으로 설정하였다. 초기 위치는 평형점이 아니며 부력이 구의 무게보다 작기 때문에 −z 방향으로 가속도가 발생해 운동이 시작된다. 자유표면 주위에는 조밀한 격자를 사용하여야 정확한 경계면을 추적할 수 있고, 계산 영역의 크기가 작아 AMR 격자 기법을 사용하지 않고 단일 격자계를 사용하였다. 계산에 사용된 전체 계산 격자수는 192×128×128이며, 균일한(uniform) 격자를 사용하였다.

[Fig. 6] Physical configuration and boundary conditions of the sphere floating on the free surface

   4.2 해석 결과

Fig. 7은 시간에 따른 구 중심의 z 방향 속도와 구 중심의 z 방향 위치를 나타낸 그림이다. 초기 구의 위치가 정확히 수면의 중앙에 위치하기 때문에 –z 방향으로 작용하는 구의 무게가 +z 방향으로 작용하는 물의 부력보다 크다. 힘의 합력이 –z 축 방향(중력방향)이기 때문에 구는 중력방향으로 가라앉기 시작한다. 구의 무게와 유체로부터 받는 힘이 평형상태에 도달하기까지 약 0.12초 동안 –z 방향으로 가속한다. 0.12초를 지난 시점부터 구의 무게보다 구가 유체로부터 받는 힘이 커지기 시작하면서 구의 속도는 감속한다. 약 0.25초에서는 구의 속도는 0이 되고 구는 가장 낮은 위치(z =−0.28m)에 도달하게 된다. 0.25초를 지난 시점부터 구는 +z 방향으로 움직이기 시작한다. 구는 상하 반복운동을 하면서 평형상태로 근접해 간다.

[Fig. 7] Results of the sphere floating on the free surface: (a) time history of the z？direction velocity, (b) time history of the center position

Fig. 8은 구가 중력방향으로 가속 운동하는 순간(t=0.63초)에서의 속도벡터와 자유표면의 모양을 𝑥z 평면에서 나타내었다. 속도벡터의 경향성을 알아보기 위해 각 방향으로 4개 중 1개의 벡터만 표시하였다. 구의 –z 방향 운동은 구 주위 유체를 밀어내게 되어 Fig. 8의 A 부근의 유체를 +z 방향으로 밀어 올리게 된다. 또한 B 부근에서는 수면의 높이가 주위보다 높아 중력에 의해 –z 방향으로 밀려 내려가게 된다. 시간이 지남에 따라 수면의 상하가 반복적으로 변동해 간다. 유동장은 구를 중심으로 대칭적으로 형성되는 것을 확인할 수 있다.

[Fig. 8] Instantaneous velocity field and free surface in xz plane during a downward movement of the sphere

Fig. 9는 시간에 따른 구가 받는 z방향 힘을 나타낸 그림이다. 구가 받는 힘은 시간이 지남에 따라 구의 무게로 수렴해 가는 것을 확인할 수 있다. 또한 구에 작용하는 힘은 구의 무게를 중심으로 진동양상을 보이는데 이로 인해 구는 상하 왕복운동을 하게된다.

[Fig. 9] Time history of the buoyant force acting on the sphere (z direction)

Fig. 10은 구가 중력방향으로 운동하는 순간의 자유표면의 모양을 나타낸 그림이다. 구의 진동에 의한 수면 wave가 방사형으로 퍼져 나가는 것을 확인할 수 있다.

[Fig. 10] Instantaneous snapshot of the free-surface during a downward movement of the sphere

5. 해양구조물 주위의 다상유동 해석

   5.1 해양구조물의 모델링

해양구조물 중에서 부유식의 일종인 반잠수식(semisubmersible) 해양구조물에 대한 다상유동 해석을 수행하였다. 본 연구에 사용된 반잠수식 해양구조물의 AutoCAD로 모델링된 형상은 Fig. 11(a)와 같다. AutoCAD 파일로 부터 모델 표면의 점(Marker particle)들을 자동적으로 추출할 수 있는 프로그램 (Choi, et al., 2011)을 이용하여 표면점들을 자동적으로 추출하였다. Fig. 11(b)는 Solver에서 인식한 형상을 Tecplot을 이용하여 나타낸 그림이며, 실제 모형과 잘 일치하는 것을 확인할 수 있다. 본 연구에 사용된 해양구조물의 크기는 Fig. 11(b)에 나타내었다. 해양구조물은 하나의 상부구조(superstructure), 2개의 하부구조(substructure), 4개의 기둥(column)으로 구성되어 있다(Fig. 11(b)).

[Fig. 11] Modeling of semisubmersible: (a) AutoCAD modeling, (b) structure of marker particles plotted by Tecplot

   5.2 계산영역 및 경계조건

전체 계산영역 및 경계조건은 Fig. 12와 같다. 전체 계산 영역의 크기는 1200m×600m×150m 이다. 입구에서는 Dirichlet 경계조건(u=2knots (=1.029m/s), 𝜐=0m/s , w =0m/s)을 사용하였고, 출구에서는 대류 경계조건을 사용하였다. 바닥면(–z방향 면)에는 점착조건(no-slip)을 사용하였고, 나머지 면에는 slip 조건을 사용하였다. 중력은 –z 방향으로 설정하였다. 해양구조물의 바닥면 위치를 좌표계 z축의 원점으로 설정하였다. 해양구조물의 바닥면과 계산영역의 바닥면 사이의 거리는 20m이고 (Fig.13), 초기 자유표면의 위치는 z =20m 이다. 계산에 사용된 물(water)과 공기(air)의 물성치는 각각 ρ_water=998.2kg/m³, ρ_air=1.225kg/m³, μ_water =1.003×10⁻³kg/(ms), μ_air =1.789× 10⁻⁵kg/(ms)이다. 해양구조물의 위치를 고정하고 입구에서 균일하게 유입되는 다상유체에 의한 수면의 변화 및 해양구조물 후류의 유동 특성을 관찰하였다.

[Fig. 12] Physical configuration and boundary conditions of the marine structure partially submerged under the free surface

[Fig. 13] Computational domain in ？z plane: blue line represents the free surface

AMR 격자구조를 사용하여 해석을 수행하였다. 4단계의 블럭을 사용하였고, 블럭들은 각 방향으로 8개의 격자로 구성되어 있다. 가장 조밀한 블럭에서의 격자 크기는 1.17m×1.17m×1.17m이다. 해양구조물 주위에는 조밀한 블럭을 배치하였고, 해양구조물에서 멀어질수록 성긴 블럭을 배치하였으며, 사용된 총 블럭수는 2958개이다(Fig. 14). 해양구조물과 다상유체 계산에 사용된 Re는 8.2×10E7 이다. 이는 입구에서의 유입속도와 물의 물성치, 해양구조물의 𝓎방향 길이로 계산하였다. 사용된 격자계가 해양구조물 주위의 경계층을 정확히 해석할 만큼 조밀하지는 못하지만, 복잡한 형상의 해양구조물과 다상유체와의 상호작용 해석을 구현할 수 있다는 점에 본 논문의 초점을 맞추었다.

[Fig. 14] Three-dimensional block distribution for the case of flow past a marine structure

   5.3 해석 결과

Fig. 15는 시간에 따라 해양구조물이 받는 힘을 나타낸 그림이다. x 방향 힘은 약 258kN을 기준으로 하여 진동하는 것을 확인할 수 있다. 𝓎 방향 힘은 양과 음의 값으로 반복적으로 진동한다. z 방향 힘은 약 396MN으로, 높은 정수압에 의하여 구조물은 +z 방향으로 힘을 받는다. 해양구조물이 받는 각 방향 힘의 평균과 rms(root mean square)를 Table 1에 나타내었다. 힘의 평균값은 z 방향이 높은 압력으로 인하여 가장 크지만, rms값은 𝓎방향이 더 크다.

[Fig. 15] Time history of forces acting on semisermercible: (a) x, ？ direction forces, (b) z direction force

[Table 1] Time mean and rms of forces

Time mean and rms of forces

Fig. 16은 시간에 따른 각 방향 힘을 FFT(Fast Fourier Transform)하여 나타낸 그림이다. 각 방향으로 해양구조물이 받는 힘이 주기적이고 않고 불규칙하기 때문에 여러 개의 주파수가 혼재되어 있다. x, 𝓎 방향의 경우 가장 지배적인 주파수는 0.00967Hz이며, 𝓎 방향의 경우 가장 지배적인 주파수는 0.01611Hz이며, 두 번째 지배적인 주파수는 0.00967Hz이다.

[Fig. 16] FFT of time history of force acting on marine structure: (a) Fx, (b) F？, (c) Fz

Fig. 17은 𝑥𝓎 평면에서의 순간적인 속도 벡터들을 나타낸 그림이다. Fig. 17의 붉은색 점원에 해당하는 부분에 대한 확대도는 Fig. 18에 나타내었다. Fig. 17~18 (a)는 하부구조물을 지나는 유동이며(z =5m), (b)는 기둥을 지나는 유동(z =20m), (c)는 상부구조물을 지나는 유동이다(z =38m). Fig 17~18 (a)와 (b)의 경우 속도 벡터들을 관찰하기 위해 상부구조물을 생략하여 나타내었다. 하부구조물을 지나는 유동에서는 다른 위치(상부구조물과 기둥)에 비해 하부구조물 주위에서 얇은 경계층이 발달하며(Fig. 18(a)), 하부구조물의 후류 끝에서 유동 박리가 발생한다(Fig. 17(a)). 기둥을 지나는 유동에서는 유동방향의 상류에 있는 기둥(Fig. 17(b)의 A, C) 에서 유동박리가 발생하며, 상류 기둥 후류에서 생성된 불규칙한 와류들이 하류 기둥(Fig. 17(b)의 B, D)에 부딪힌다. 하류 기둥을 지나는 유체들은 후류에서 상호작용으로 유체들이 섞이는 것을 확인할 수 있으며(Fig. 17(b)의 E) 일부 유체들이 상류와 하류 기둥 사이에서 정체되어 있는 것을 확인할 수 있다(Fig. 18(b)). 상부구조물을 지나는 공기 유동의 경우 상부구조물의 상류 쪽 모서리에서 유동 박리가 발생하며(Fig. 17(c)), 상부구조물의 옆면에서 와류들을 관찰할 수 있다(Fig. 18(c)).

[Fig. 17] Instantaneous velocity plot in ？？ plane: (a) z=5m, (b) z=20m, (c) z=38m

[Fig. 18] Enlarged instantaneous velocity plot in ？？ plane: (a) z=5m, (b) z=20m, (c) z=38m

Fig. 19는 순간적인 해양구조물 주위의 수면의 모양을 나타낸 그림이다. 물과 공기의 부피비를 나타내는 용적함수 f가 0.5인 등위면(iso-surface)을 이용하여 수면의 모양을 형상화하였다. 입구에서 유입된 평평한(flat) 상태의 수면은 상류 기둥을 만날 때까지 평평한 상태를 유지하다가 상류 기둥에서 유동 박리가 발생하는 위치 부근부터 변화를 보이기 시작한다. 수면의 진폭은 위치마다 다르지만, 약 +3.2m ∼−2m 사이에서 중력방향으로 진동하였다. 상류 기둥 사이 후류(Fig. 19의 A)에서는 수면의 진동이 거의 발생하지 않고 평평한 상태를 유지하지만, 하류 기둥 사이 후류(Fig. 19의 B)에서는 +𝓎 방향 구조물과 –𝓎 방향 구조물 후류 사이의 상호작용으로 인하여 복잡한 형태의 수면 진동이 발생한다. 유체가 유동 방향으로 진행하면서 유동 불안정성이 증가하여 하류 기둥 후류를 지나면서 서로 섞이게 되는데(Fig. 17(b)) 이로 인하여 하류 기둥 사이 후류에서 수면의 진동이 발생한다. 상류 기둥 후류에서 생성된 복잡한 형태의 수면은 주유동 방향으로 흘러가면서 하류 기둥에 부딪히게 된다. 상대적으로 상류 기둥과 하류 기둥 사이에서 발생하는 수면의 진동 진폭(Fig. 20의 C)이 하류 기둥 후류(Fig. 20의 D)에서 발생하는 수면의 진동 진폭보다 더 크다. 상류 기둥과 하류 기둥 사이에서 일부의 유체가 기둥들에 의해 갇히게 되는데, 이로 인해 z 방향으로 상승, 하강하는 유체가 발생하여 수면의 진동 진폭이 커지게 된다.

[Fig. 19] Instantaneous snapshot of the free-surface around the marine structure

[Fig. 20] Instantaneous snapshot of the free-surface around the marine structure

Fig. 21은 𝑥z 평면에서의 순간적은 속도 벡터들을 나타낸 그림이다. Fig. 21(a)는 𝓎=0인 평면에서의 속도 벡터들을 나타내었다. 상부구조물의 –z 방향의 면을 지나는 유체는 얇은 경계층을 이루면서 빠져나가지만 +z 방향의 면을 지나는 유체는 상부구조물의 +z 방향 면에서 재순환구역이 발생한다. 상부구조물을 완전히 지나서부터는 +z 방향 면을 지나는 유체와 –z 방향 면을 지나는 유체가 상부구조물 후류에서 섞인다. Fig. 21(b)는 𝓎=−30m인 평면에서의 속도벡터들이다. 하부구조물의 윗부분에서 유동 박리가 발생하며 상류 쪽 기둥사이에서 재순환구역이 발생한다. 상류, 하류 기둥 사이에서 유체들이 정체되어 있는 것을 관찰 할 수 있다.

[Fig. 21] Instantaneous velocity plot in ？z plane: (a) ？=0m, (b) ？=？30m

6. 결 론

해양구조물과 같이 복잡한 형상을 갖는 구조물 주위의 유동 해석에 있어, 구조물 주위에 조밀한 격자의 사용은 수치해석의 정확성을 좌우하는 중요한 요소이다. 그러나 조밀한 격자의 사용으로 인한 수치해석 비용의 급격한 증가는 특히 물체와 유체 사이에 상호작용을 고려해야 하는 VIV 해석에 있어서 최대의 난제중 하나이다. 물체 주위 유동 해석분야에서 기존에 주로 사용된 비정렬 격자계 혹은 Body-fitted 좌표계 기반의 해석은 격자 문제로 인한 한계가 존재한다. 본 연구에서는 이러한 격자 문제를 해결하기 위해 직교 좌표계 기반의 AMR 격자구조를 사용하였다. 구역에서는 성긴 격자를 사용하고, 물리량이 급변하는 곳에서는 조밀한 격자를 사용함으로써 효율적이고 정확한 계산을 수행할 수 있었다. AMR 격자구조에서 출발하여 기존 VIV 해석에 사용되고 있는 IBM 및 높은 레이놀즈 수의 난류를 정확히 해석하기 위한 LES 해석 기술을 적용하였고, VOF 기법을 적용하여 자유표면을 포함하는 다상 유체와 구조물과의 상호작용을 정확히 해석할 수 있는 코드를 개발하였다. 또한 AutoCAD로 모델링된 파일에서 표면점들을 추출하여 복잡한 형상을 쉽게 구현할 수 있는 기능을 추가하였다. 기존의 IBM에서는 복잡한 형상의 물체를 구현하기 위해서 물체 표면을 나타내는 방정식을 이용하였지만, 이런 방법의 경우에는 간단한 형상에만 적용할 수 있으며, 형상이 복잡해지면 정확한 물체의 표면을 구현 하지 못한다는 어려움이 있다. 이러한 문제를 해결하기 위해서 본 연구에서는 임의의 형상을 쉽게 구현하기 위해 표면점(marker particle)을 이용하는 방법을 채택하였다.

개발된 코드를 이용하여 단상 유동 및 다상 유동에 대한 검증을 수행하여 코드의 타당성을 확보하였다. 완성된 코드를 이용하여 자유표면에 떠 있는 구의 동적 거동과 수면의 변화에 대한 연구와 해양구조물 주위 수면의 변화 및 구조물이 받는 힘에 대한 연구를 수행하였다.

참고문헌

1. Choi C.B., Jang Y.J., Lee K., Yang K.S. 2011 Fire dynamic simulator with marker particles [Journal of Computational Fluids Engineering] Vol.16 P.55-63
2. Hall K.C. 1994 Eigenanalysis of unsteady flows about airfoils, cascades, and wings [AIAA Journal] Vol.32 P.2426-2432
3. Hirt C.W., Nichols B.D. 1981 Volume of Fluid(VOF) method for the dynamics of free boundaries [Journal of Computational Physics] Vol.39 P.201-225
4. Jadic I., So R.M.C., Mignolet M.P. 1998 Analysis of fluid-structure interactions using a time marching technique [Journal of Fluids and Structures] Vol.12 P.631-654
5. Kim J., Moin P. 1985 Application of a fractional-step method to incompressible Navier-Stokes equations [Journal of Computational Physics] Vol.59 P.308-323
6. Koshizuka S., Oka Y. 1996 Moving-particle semi-implicit method for fragmentation of incompressible fluid [Nuclear science and engineering] Vol.123 P.421-434
7. Martin J.C., Moyce W.J. 1952 An experimental study of the collapse of liquid columns on a rigid horizontal plane [Philosophical transactions of the royal society A] Vol.244 P.312-324
8. Meneveau C., Lund T., Cabot W. 1996 A lagrangian dynamic subgrid-scale model of turbulence [Journal of fluid mechanics] Vol.319 P.353-385
9. Park J.C, Kim M.H, Miyata H., Chun H.H. 2003 Fully nonlinear numerical wave tank(NWT) simulations and wave run-up prediction around 3-D structures [Ocean Engineering] Vol.30 P.1969-1996
10. Pilliod J.E. 1992 An analysis of piecewise linear interface reconstruction algorithms for volume-of-flud methods. M. S. Thesis
11. Pilliod J.E., Puckett E.G. 2004 Second-order accurate volume-of-fluid algorithms for tracking material interfaces [Journal of Computational Physics] Vol.199 P.465-502
12. Robertson I., Li L., Sherwin S.J., Bearman P.W. 2003 A numerical study of rotational and transverse galloping rectangular bodies [Journal of Fluids and Structures] Vol.17 P.681-699
13. Vanella M., Rabenold P., Balaras E. 2010 A direct-forcing embedded-boundary method with adaptive mesh refinement for fluid-structure interaction problems [Journal of Computational Physics] Vol.229 P.6427-6449
14. Yang J., Balaras E. 2006 An embedded-boundary formulation for large-eddy simulation of turbulent flows interacting with moving boundaries [Journal of Computational Physics] Vol.215 P.12-40

이미지 / 테이블

[ ]
[ ]
[ ]
[ ]
[ Fig. 1 ] AMR block system
[ ]
[ ]
[ ]
[ Fig. 2 ] Physical configuration and computational domain of the free-falling sphere problem
[ ]
[ Fig. 3 ] Results of the free-falling sphere problem: (a) time history of the velocity, (b) time history of the drag coefficient
[ Fig. 4 ] Physical configuration of a broken water-dam problem
[ Fig. 5 ] Comparison of length of the surge front(z/？) against time
[ Fig. 6 ] Physical configuration and boundary conditions of the sphere floating on the free surface
[ Fig. 7 ] Results of the sphere floating on the free surface: (a) time history of the z？direction velocity, (b) time history of the center position
[ Fig. 8 ] Instantaneous velocity field and free surface in xz plane during a downward movement of the sphere
[ Fig. 9 ] Time history of the buoyant force acting on the sphere (z direction)
[ Fig. 10 ] Instantaneous snapshot of the free-surface during a downward movement of the sphere
[ Fig. 11 ] Modeling of semisubmersible: (a) AutoCAD modeling, (b) structure of marker particles plotted by Tecplot
[ Fig. 12 ] Physical configuration and boundary conditions of the marine structure partially submerged under the free surface
[ Fig. 13 ] Computational domain in ？z plane: blue line represents the free surface
[ Fig. 14 ] Three-dimensional block distribution for the case of flow past a marine structure
[ Fig. 15 ] Time history of forces acting on semisermercible: (a) x, ？ direction forces, (b) z direction force
[ Table 1 ] Time mean and rms of forces
[ Fig. 16 ] FFT of time history of force acting on marine structure: (a) Fx, (b) F？, (c) Fz
[ Fig. 17 ] Instantaneous velocity plot in ？？ plane: (a) z=5m, (b) z=20m, (c) z=38m
[ Fig. 18 ] Enlarged instantaneous velocity plot in ？？ plane: (a) z=5m, (b) z=20m, (c) z=38m
[ Fig. 19 ] Instantaneous snapshot of the free-surface around the marine structure
[ Fig. 20 ] Instantaneous snapshot of the free-surface around the marine structure
[ Fig. 21 ] Instantaneous velocity plot in ？z plane: (a) ？=0m, (b) ？=？30m