2.1 지배방정식

본 연구에 적용된 오픈소스 CFD 코드인 REEF3D(Bihs et al., 2016)는 비압축성 혼상유동(Two-Phase Flow)을 해석하기 위해 3차원 RANS(Reynolds-Averaged Navier-Stokes equations) 방정식이 운동량 및 질량보존을 규정하기 위한 연속방정식과 함께 적용된다. 지배방정식은 다음의 식(1)과 (2)에 나타내는 연속방정식과 운동량보존방정식으로 구성된다.

여기서, u는 유속벡터, ρ는 밀도, p는 압력, ν는 동점성 및 난류와점성계수, g는 중력가속도이다. 난류모델은 Wilcox(1994)에 의해 제안된 난류운동에너지 k(Turbulent kinetic energy)와 난류비소산율 ω(Specific dissipation rate)로 구성되는 2-방정식(Two-equations)의 k-ω 모델이 적용되어, 확산항에 동점성 및 난류와점성을 고려하여 계산된다.

물과 공기에 대한 혼상유동(Two-Phase Flow)의 경계면인 자유수면의 해석에는 Level-Set Method(Osher and Sethian, 1988)가 적용되며, 각 상(Phase)의 경계면 위치는 부호가 부여된 평활화된 거리함수 Φ에 의해 표현된다. Φ는 경계면에 가장 가까운 거리를 제공하고, 부호의 변화에 따라 각각의 유체에 대한 상(Phase)을 구별한다. 즉, 물과 공기의 경계면을 0으로 하여 물의 영역은 경계면을 기준으로 양(+)의 부호, 공기의 영역은 음(-)의 부호를 거리에 따라 부여한다. 거리함수 Φ 및 경계면에 대한 시간발전의 해석에는 식(3)의 이류방정식이 적용된다.

거리함수 Φ가 흐름에 의해 이류되면 부호 거리 속성이 상실됨에 따라 이 속성을 유지하고 질량 보존을 보장하기 위해 매시간 단계마다 다시 초기화되어야 한다. 따라서, Sussman et al.(1994)의 PDE(Partial Differential Equation) 기반 재초기화 방정식이 적용된다. 각 유체에 대한 밀도(ρ)와 점성(ν)은 식(4)와 식(5)를 사용하여 계산된다.

여기서, H는 Heaviside 함수, 아래첨자 w, a는 각각 물과 공기를 나타낸다.

경계면에서의 밀도와 점성의 급격한 변화는 수치적 불안정을 야기할 수 있다. 따라서 경계면에는 격자크기에 비례하는 경계면 두께 ε = 2.1△x을 부여하고, 해당 영역의 평활화를 위한 Heaviside 함수 H(식(6))를 적용한다.

식(1)~(3)의 지배방정식은 직교 교호격자를 적용한 유한차분법(Finite Difference Method, FDM)에 의해 이산방정식이 구성되며, 공간이산화에는 5차 WENO Scheme(Jiang and Peng, 2000; Jiang and Shu, 1996)이 적용되고 확산항은 2차 중앙차분으로 이산화된다. 시간이산화에는 3차 TVD Runge-Kutta Scheme(Shu and Osher, 1988)이 적용되고, 난류모델은 Implicit Euler method가 적용되어 CFL 조건에서 격자 크기에 대한 의존성을 효과적으로 제거한다. 압력은 Incremental pressure-correction algorithm(Timmermans et al., 1996)을 식(1)~(2)에 적용하여 계산된다. 이때 3차 TVD Runge-Kutta Scheme의 각 k번째 단계별 유속장 u(*)는 식(7)을 이용하여 이전 단계의 압력구배를 사용하여 예측된다.

여기서, α_k = 1.0, 1/4, 2/3, β_k = 0.0, 3/4, 1/3 그리고 k = 1, 2, 3을 나타낸다. 이후, HYPRE solver library(Center for Applied Scientific Computing, 2006)의 PF_MG multi-grid solver(van der Vorst, 1992)로 사전조정과 함께 Fully parallelized BiCGStab algorithm을 사용하여 압력 보정항 p_corr에 대해 식(8)의 Poisson 압력방정식이 계산된다.

그리고 압력장 및 유속장은 최종적으로 식(9)와 식(10)을 이용하여 계산된다.

2.2 유체와 구조물의 상호작용 해석(Fluid-Structure Interaction, FSI)

CFD 모델에서 운동하는 구조물과 상호작용하는 이상유동(Two-Phase Flow)을 수치해석하는 FSI의 기법은 격자체계의 관점에서 크게 Moving grid method와 Fixed grid method로 분류할 수 있다. Moving grid method의 경우 구조물의 형상에 따라 계산격자를 재구성함으로 구조물 인근에서 정확한 결과를 얻을 수 있지만, 중첩 격자영역의 보간 및 격자 재구성에 따른 계산비용의 증가 등의 문제점도 내포하고 있다. 반면에 Fixed grid method의 대표적인 기법으로 Peskin(1972)가 제안한 Immersed Boundary Method는 배경격자가 변형되지 않고 복잡한 형상의 구조물 이동 및 재배치시에 격자를 재생성할 필요가 없는 장점이 있다. 이를 기반으로 Fadlun et al.(2000), Uhlmann(2005), Yang and Stern(2015), Yang (2018) 등의 연구에 의해 Direct forcing method/Immersed Boundary Method의 개발 및 추가적인 발전이 수행되었다.

본 연구에 적용된 오픈소스 CFD는 유체중에 운동하는 강성 구조물의 해석방법으로 구조물의 내부 영역을 주위와 동일한 유체로 채우고 구조물의 운동을 위하여 내부의 유체 영역에 가상의 외력을 작용시키는 Yang(2018)에 의해 제안된 Continuous direct forcing Method/Immersed Boundary Method가 적용되었다(Martin at al. 2021). 유체영역을 점유하는 구조물은 Eulerian 격자영역에서 Lagrangian point를 통해 구조물 형상을 표현하게 되며, 자유수면의 해석과 동일한 Level set method가 적용된다. 구조물은 삼각형 표면격자(Triangular surface mesh)로 구성된 STL(STereoLitography format)형식으로 솔버에 입력되고, 입력된 STL 정보는 Ray-tracing algorithm을 적용하여 구조물의 내부-외부 정보를 Level set 함수 Φ_s로 표현한다. 생성된 Level set 함수 Φ_s는 식(4)~(5)를 확장하여 유체 및 구조물 영역을 구분하는데 적용된다.

여기서, 하첨자 s는 구조물을 나타낸다.

Yang(2018)에 의해 제안된 연속방정식과 운동량보존방정식은 식(13)과 식(14)로 각각 구성된다.

여기서, 외력항 f는 식(15)로 주어진다.

여기서, P(u)는 유체의 유속장을 Lagrangian point로 이루어진 구조물 내부로 투영하는 유속장을 나타낸다.

식(14)를 식(1)~(2)와 비교하면, 유일한 차이는 외력항 f와 확산항임을 알 수 있다. 따라서, 유체에서 구조물로의 전이를 나타내는 Heaviside 함수 H(Φ_s)와 구조물 경계면에서의 외력항 f를 적용하여 전체 영역에서 단일 방정식을 적용하여 계산된다. 이산화 단계에서 외력항 f는 새로운 시간단계 n + 1에서 식(16)로 나타낼 수 있다.

시간단계 n + 1에서의 P(u⁽ⁿ⁾)은 미지수이기 때문에 유효한 근사값 P(u⁽ⁿ⁾) = u⁽ⁿ⁾)를 적용하고 이전 시간 단계의 압력을 근사치로 취하면 식(17)로 나타낼 수 있다.

그런 다음 식(7)을 적용하여 식(18)로 나타낼 수 있다.

식(18)의 계산에서 식(7)을 외력항 f를 고려하지 않고 먼저 실행하여 P(u⁽ⁿ⁺¹⁾)의 근사치로 P(u^(*))를 취하면 f^(*)는 식(19)로 산정할 수 있다.

그리고, f^(*)는 식(8)의 Poisson 압력방정식을 풀기 전에 식(7)에 추가하여 계산한다. 여기서, 구조물 내부로 투영되는 유속장 P(u^(*))는 Lagrangian point로 이루어진 구조물의 운동이 병진운동과 회전운동의 조합으로 구성된다는 가정에 기초하여 해석을 수행하여 산정된다.

구조물의 병진운동은 식(20)과 같이 뉴턴의 제2법칙으로 표현된다.

여기서, F_X는 구조물의 표면력, ρ_s는 구조물의 밀도, V_s는 구조물의 부피, X˙s는 구조물 무게중심위치에서의 속도이다.

구조물의 회전운동은 고정좌표계에서 오일러 각도(Tait-Bryan angles) roll(Φ), pitch(Θ), yaw(Ψ)와 관련된 4개의 오일러 매개변수(quaternion) e = (e₀, e₁, e₂, e₃)^T를 적용하여 해석된다.

오일러 각도(Tait-Bryan angles)는 다음의 관계식에 의해 구할 수 있다.

여기서, 4개의 오일러 매개변수(quaternion)는 e^Te = 1의 조건으로 결정된다. 그리고, Shivarama and Fahren-thold(2004)이 제안한 Hamiltonian formulation을 적용하여 구조물 회전운동의 모멘텀 벡터를 산정한다. 이후, 구조물의 표면력(F_X)과 모멘트(M_X)는 구조물의 모든 삼각형 표면격자(Triangular surface mesh)에서의 유체력을 적분하여 계산된다.

여기서, n은 표면법선벡터, τ는 점성응력텐서, N은 삼각형 표면격자의 개수, r은 거리벡터를 나타낸다. 그리고, 힌지를 중심으로 회전하는 진자운동의 구현을 위해 Spring 계류 조건을 적용하였으며, 표면력(F_X)과 모멘트(M_X)에 각각 계류력 F_M과 무게중심에서 힌지위치까지의 거리벡터가 고려된 모멘트 r_CM × F_M를 부가함으로써 구현된다. 계산된 모멘트로부터 구조물 내부로 투영되는 유속장 P(u^(*))은 식(25)로부터 산정할 수 있다.

여기서, X˙s는 구조물 무게중심위치에서의 속도이고 ω_s는 구조물의 회전운동속도이다.

2.3 파랑 조파 및 흡수기법

파랑의 조파기법에는 Relaxation method(Jacobsen et al., 2012), 흡수기법에는 Active wave absorption method(Schäffer and Klopman, 2000)를 적용하였다.

수치파동수조의 조파경계에 적용된 Relaxation method는 파동이 생성되는 입구에 식(26)의 이완영역(Relaxation zone)을 적용하여 활성화된다.

여기서, x~는 이완영역의 거리에 맞게 조정되는 무차원거리를 나타내고, u는 수평방향 수립자유속, w는 수직방향 수립자유속, ϕ는 Level set 함수, p는 x~위치의 압력, Γ(x~ )는 Jacobsen et al.(2012)에 의해 제안된 파동제어를 위한 이완함수(Relaxation function)이며 식(26)과 같다.

조파경계에서 파동의 조파는 이완함수 Γ(x~)를 사용하여 파동이론으로부터 계산된 유속, 자유수면 및 압력에 대한 이론해에서 수치해로의 점진적인 변화를 통해 도입되며, 조파경계로의 재반사파는 수치해에서 이론해로의 변환을 통해 흡수된다. 이러한 방법을 통해 조파경계로 되돌아오는 반사파는 흡수되어 조파경계에서 발생하는 간섭을 최소화한다.

수치파동수조의 흡수경계 처리는 반사파를 상쇄하기 위해 흡수되는 파동유속의 반대 부호로 수정된 유속을 수치파동수조 출구에서 규정하는 Active wave absorption method를 적용하였다. 2차원 파동의 Airy 파동이론에서 천수조건을 가정하면 수평 수립자 유속은 수직축을 따라 일정하며, Active wave absorption method 경계조건은 식(27)과 같이 구현된다.

여기서, U(t)는 수정된 유속, η(t)는 흡수경계의 자유수면변위, d는 수심을 나타낸다.

4.1 수치파동수조 구성

Fig. 3에 나타내는 바와 같이 수치해석을 위한 수치파동수조의 구성은 수리모형실험 결과를 정량적으로 비교하기 위해 구조물 및 입사파랑 제원은 수리모형실험과 동일하게 설정하였다. 단, 계산의 효율성과 적용된 수치모델이 무반사 조파 기능이 포함되어 있음을 고려하여 수치파동수조의 길이는 13.96 m, 폭은 0.05 m로 감소시켜 설정하며, 높이는 1.0 m의 영역으로 설정하였다. 격자크기는 △x = 0.005~0.03 m, △y = 0.01m, △z = 0.01m의 가변격자체계로 구성하여, 적용된 격자의 총 개수는 423,594개이다. 구성된 격자에 대한 Courant 상수는 0.3 이하가 되도록 자동으로 시간 간격을 조정하며 수치해석이 수행되고 계산의 효율을 위해 영역분할에 의한 MPI 병렬계산을 적용하였다. 수조의 바닥면 및 구조물의 표면경계는 No-slip 조건을 적용하였다. 대상파랑에 대한 수치조파에는 Cnoidal파 이론을 적용하였으며, 조파(1L)경계는 Relaxation Method, 흡수경계는 Active wave absorption method(AWA)를 적용하였다.

Fig. 3.

Definition sketch of numerical wave tank.

4.2 역삼각형 플랩 구성

역삼각형 플랩은 상부 길이 15 cm, 폭 5 cm, 높이 40 cm 형상의 STL 형식으로 구성되어 솔버에 입력되고, 입력된 Fig. 4(a)의 STL 정보는 Ray-tracing algorithm을 적용하여 구조물의 내부-외부 정보를 Fig. 4(b)의 Level set 함수 Φ_s로 표현하여 유체 및 구조물 영역을 구분하게 된다. CFD 모델 내에서 Φ_s = 1인 경우를 구조물로 취급하게 된다. 힌지조건으로는 Spring 계류를 적용하였으며, 힌지위치에서 역삼각형 플랩의 하단까지 Spring 계류라인으로 연결하여 힌지를 중심으로 회전하는 진자운동을 구현하였다. 따라서, 역삼각형 플랩의 거동은 2자유도 운동(surge, pitch)으로 제한된다.

Fig. 4.

Liquid and solid level-set function distributions around the OWSC model.

4.3 자유수면변위

Fig. 5~6은 Table 1의 계산조건에 대하여 역삼각형 플랩의 전면과 배후의 자유수면변위에 대해 수리실험 및 수치해석 결과를 비교하여 나타낸 것이다. 자유수면변위의 측정위치는 Fig. 3에 제시한 바와 같이 역삼각형 플랩의 전면 3점(W1(x = 9.00 m), W2(x = 9.50 m) 및 W3(x = 9.70 m))과 후면 1점(W4(x = 12.50 m))에 위치한다.

Fig. 5.

Comparison of water surface elevations between simulated and measured (Case1~4).

Fig. 6.

Comparison of water surface elevations between simulated and measured (Case5~8).

Fig. 5(a)와 Fig. 5(b)에 나타내는 주기가 비교적 작은 Case1(T = 1.6 s) 및 Case2(T = 2.0 s)의 수면변동은 역삼각형 플랩의 전면(W1, W2, W3)과 후면(W4)에서 모두 비교적 선형적인 파형을 유지하는 것을 보인다. 이는 진행하는 파동과 역삼각형 플랩의 상호작용에 의한 거동이 파동장을 크게 교란시키지 않고 연행되고 있음을 보여준다.

다음으로 Fig. 5(c)에 나타내는 주기가 소폭 증가하는 Case3(T = 2.4 s)의 수면변동에서는 파봉에 비해 파곡 부분이 평탄해지며 경사가 커지는 파형의 상하 비대칭성이 나타나는 것이 관찰된다. 그리고 점차 주기가 증가하는 Fig. 5(d) 및 Fig. 6(a)~Fig. 6(d)에 보이는 Case4(T = 2.8 s)~Case8(T = 4.0 s)의 수면변동은 파형의 상하·전후 비대칭성 및 비선형성이 크게 나타나고 있고, 특히 파형 사이에 고주파수 파형성분이 발달됨을 확인할 수 있다. 이는 역삼각형 플랩의 부력에 의한 복원력으로 인해 발생되는 회전운동 및 관성운동이 진행하는 파동과 상호작용하여 생성되는 반사파동의 영향으로 판단된다. 전반적으로 주기가 길어질수록 비대칭성과 비선형성이 강하고 고주파수 파형성분이 발달되는 수위변동이 형성되는 것을 알 수 있다.

전반적으로 주기가 길어질수록 비대칭성과 비선형성이 강하고 고주파수 파형성분이 발달되는 수위변동이 형성되는 것을 알 수 있으며, 수리실험과 수치해석결과의 대응성은 자유수면변위의 진폭과 위상에서 약간의 편차를 보이고 있으나, 시간발전에 따른 역삼각형 플랩 주변에서 발생하는 일련의 파랑변형과정을 양호하게 재현하고 있다는 것을 확인할 수 있다.

4.4 플랩회전각도

Fig. 7은 역삼각형 플랩의 x 축 방향 회전운동(pitch)에 의한 회전각도(Flap Rotation angle)에 대해 수리실험 및 수치해석 결과를 비교하여 나타낸 것이다. 그래프의 세로축에 나타내는 회전각도에서 양의 값(+)은 시계방향 회전, 음의 값(-)은 반시계방향 회전을 나타낸다. 전반적으로 회전각도의 시간발전에 따른 양상은 모든 Case에서 비교적 선형적으로 나타나고, Case1(T = 1.6 s)에서 Case8(T = 4.0 s)로 주기가 증가함에 따라 시계방향 및 반시계방향의 회전각도는 증가되어 나타난다.

Fig. 7.

Comparison of flap rotation angle between simulated and measured.

Fig. 8은 시계방향 및 반시계방향 회전각도의 최대치를 나타낸 것으로 시계방향 회전의 경우, Case1(T = 1.6 s)과 Case8(T = 4.0 s)에서 수리실험 결과는 13°와 29°, 수치해석 결과는 16°와 27°로 각각 16°, 11°의 회전각도 증가가 나타난다. 그리고 반시계방향 회전의 경우, 수리실험 결과는 -13°와 -29°, 수치해석 결과는 -10°와 -22°로 각각 16°, 12°의 회전각도 증가가 나타난다. 이때 회전각도는 Case5(T = 3.0 s)까지 지속적으로 증가하다가 Case6(T = 3.2 s)을 기점으로 증가폭이 둔화되어 Case8(T = 4.0 s)로 갈수록 일정각도로 수렴되어 나타나는 결과를 보인다. 이는 본 연구에서 적용하는 실험 Case가 주기의 변화만을 고려하므로, 실험 Case의 변화로 인한 수평유속의 증가는 크지 않고 역삼각형 플랩에 작용하는 외력의 지속시간 증가와 복원력과의 평형으로 야기되는 결과로 판단된다.

Fig. 8.

Comparison of maximum flap rotation angle between simulated and measured.

Fig. 7과 Fig. 8에서 보이는 바와 같이, 수리실험과 수치해석결과의 대응성은 Case1(T = 1.6 s)에서 Case8(T = 4.0 s)로 주기가 증가함에 따라 반시계방향 회전각도에서 다소 차이를 나타내지만, 시간발전에 따른 회전각도의 진폭과 위상의 재현성은 양호한 것으로 판단된다.

4.5 자유수면변위 및 플랩회전각도의 오차평가

자유수면변위 및 플랩회전각도에 대한 수치해석결과와 수리실험결과의 정량적인 오차평가를 위해 정규화된 평균 제곱근 오차(Normalized Root Mean Squared Error, 이하 NRMSE)를 비교하였다. NRMSE의 계산은 식(28)을 적용하였고, Case별 NRMSE를 Table 2에 나타내었다.

Table 2.

Results of NRMSE for water surface elevations and flap rotation angle

여기서, S_i는 i번째 수치해석결과, O_i는 i번째 수리실험결과, n은 총 데이터수이다.

Table 2에 보이는 바와 같이 수리실험과 수치해석결과의 대응성은 자유수면변위의 경우 4.00~16.96%, 플랩회전각도의 경우 3.67~15.70%의 NRMSE 수준에서 재현성을 확보하는 것을 확인할 수 있다.

4.6 플랩거동특성

본 수치해석의 결과 중 Case5(T = 3.0 s)의 한 주기의 파동(t/Ti = 1/10)에 대한 역삼각형 플랩 주변의 유속벡터 및 수면 변동에 대한 Snapshot을 Fig. 9~10에 제시한다. 이때, 입사파랑은 좌측에서 우측으로 진행되며, 그림에서 3차원 Snapshot은 역삼각형 플랩의 회전운동에 대한 거동특성을 명확하게 파악하기 위해 y 축을 10배 확장시켜 표현한 것이다.

Fig. 9.

Snapshot of the velocity vectors and free surface elevation around the OWSC for Case5 (t/T = 1/10~5/10).

Fig. 10.

Snapshot of the velocity vectors and free surface elevation around the OWSC for Case5 (t/T = 6/10~10/10).

Fig. 9에서 내습하는 파랑의 파봉이 역삼각형 플랩으로 접근함에 따라 유속벡터로부터 확인되는 우측방향으로 작용하는 수평유속에 의해 역삼각형 플랩은 시계방향 회전운동을 시작하는 것을 확인할 수 있다(Fig. 9(a), (b) 참조). 파랑이 계속 진행하여 역삼각형 플랩이 파곡과 파봉의 중간지점에 위치하였을때 시계방향 회전각도가 최대로 나타나는 것이 Fig. 9(c)를 통해 확인된다. 그리고 Fig. 9(d), (e)에서 작용하는 수평유속이 좌측방향으로 전환되며 부력에 의한 복원력으로 인해 원위치로 복귀하려는 반시계방향 회전운동이 이행되는 과정을 확인할 수 있다.

Fig. 10에서 역삼각형 플랩이 파곡에 위치하는 순간에 회전각도가 원위치로 되는 것이 Fig. 10(a)을 통해 확인되고, 관성력과 유속벡터로부터 확인되는 좌측방향의 수평유속에 의해 역삼각형 플랩의 반시계방향 회전운동이 가속된다(Fig. 10(b), (c) 참조). 시계방향 회전운동과 마찬가지로 역삼각형 플랩이 파봉과 파곡의 중간지점에 위치하였을때 반시계방향 회전각도가 최대로 나타나는 것이 Fig. 10(c)를 통해 확인된다. 그리고 Fig. 10(d), (e)에서 작용하는 수평유속이 우측방향으로 전환되며 부력에 의한 복원력으로 인해 원위치로 복귀하려는 시계방향 회전운동이 이행되는 과정을 확인할 수 있다. 이상의 결과로부터 파랑과 역삼각형 플랩의 상호작용에 의한 유체역학적 거동 특성을 확인할 수 있다.