비특이화 간접경계적분방정식방법을 이용한 2차원 수치수조 개발 및 적용

Development and Application of Two-Dimensional Numerical Tank using Desingularized Indirect Boundary Integral Equation Method

Article information

J. Ocean Eng. Technol. 2018;32(6):447-457
Publication date (electronic) : 2018 December 27
doi : https://doi.org/10.26748/KSOE.2018.32.6.447
*Korea Research Institute of Ships and Ocean Engineering, Daejon, Korea
오승훈,*orcid_icon, 조석규*, 정동호*, 성홍근*
*한국해양과학기술원 부설 선박해양플랜트연구소
Corresponding author Seunghoon, Oh: +82-42-866-3921, carot541@kriso.re.kr
It is noted that this paper is revised edition based on proceedings of KAOST 2018 in Jeju.
Received 2018 August 6; Revised 2018 November 16; Accepted 2018 November 19.

Trans Abstract

In this study, a two-dimensional fully nonlinear transient wave numerical tank was developed using a desingularized indirect boundary integral equation method. The desingularized indirect boundary integral equation method is simpler and faster than the conventional boundary element method because special treatment is not required to compute the boundary integral. Numerical simulations were carried out in the time domain using the fourth order Runge-Kutta method. A mixed Eulerian-Lagrangian approach was adapted to reconstruct the free surface at each time step. A numerical damping zone was used to minimize the reflective wave in the downstream region. The interpolating method of a Gaussian radial basis function-type artificial neural network was used to calculate the gradient of the free surface elevation without element connectivity. The desingularized indirect boundary integral equation using an isolated point source and radial basis function has no need for information about the element connectivity and is a meshless method that is numerically more flexible. In order to validate the accuracy of the numerical wave tank based on the desingularized indirect boundary integral equation method and meshless technique, several numerical simulations were carried out. First, a comparison with numerical results according to the type of desingularized source was carried out and confirmed that continuous line sources can be replaced by simply isolated sources. In addition, a propagation simulation of a 2nd-order Stokes wave was carried out and compared with an analytical solution. Finally, simulations of propagating waves in shallow water and propagating waves over a submerged bar were also carried and compared with published data.

1. 서 론

선박 및 해양구조물의 안전한 설계, 건조 및 운용은 해양공학 분야의 주요 관심사이다. 이는 선박 및 해양구조물에 승선한 승무원 및 선적된 화물의 안전과 구조물의 안정성 확보에 직결되기 때문이다. 특히 극한 해양환경에서 선박 및 해양구조물의 거동은 설계 및 운용성 기준을 결정하는 주요한 인자가 된다. 그 중 파랑은 구조물의 동적 거동에 가장 큰 영향을 주는 환경요소이다.

따라서 해상의 파랑현상을 정확히 이해하고 설계파를 산정하는 것이 필수적이다. 파도의 생성 및 전파는 파랑현상에서 이해해야할 중요한 요소 중 하나이다. 특히 해안의 수심변화에 따라 파의 진폭, 파장 및 파향 등과 같은 파의 주요특성들이 변화하며 파정의 경사가 급하고 파저가 완만한 비대칭 파형으로 대표되는 비선형성이 나타난다. 이와 같이 시간에 따라 변화하는 비선형적인 파의 생성 및 전파 현상을 잘 이해하기 위한 수조시험 또는 수치수조의 역할이 증대되고 있다. 복잡한 불규칙 해상을 시험적 또는 수치적으로 구현하는 것은 쉽지 않다. 수십 년 전부터 조파수조와 같은 시험 설비들이 제작되어 파의 특성은 물론 이에 기진된 힘과 구조물의 운동 등에 대한 분석이 진행되어 현상에 대한 물리적인 이해를 돕고 있다. 현재 해양과학기술원 부설 선박해양플랜트연구소(Korea Research Institute of Ships and Ocean Engineering, KRISO)에서는 해양공학수조(Ocean engineering basin, OEB)를 운용하고 있으며, 체계적인 심해 및 천해 조파시험이 가능한 심해공학수조(Deepwater ocean engineering basin, DOEB; 수심 15m; 수심 조절 가능)를 건설 중에 있다. 조파시험의 경우 시간과 비용이 많이 소요되므로 수조 일정이 제한적 이다. 그러므로 조파수조의 시험조건을 선정하고 시험결과를 해석하는 효율적인 수단으로 수치수조가 활용될 수 있다. 극한조건의 파랑은 비선형성이 강하기 때문에 파 간의 상호작용을 구현하는데 선형이론 또는 2차 비선형 이론으로는 한계가 있는 것으로 알려져 있다. Ducrozet et al.(2012)는 2차 비선형 수치수조, 고차 스펙트럴 수치수조 및 시험수조에서 불규칙 파를 구현하여 비교하였으며 2차 비선형 수치수조는 수조시험의 파고를 정확히 구현하지 못하는 것으로 확인하였다. 앞서 언급한바와 같이 심해에서 천해로 천이하는 과정에서 파의 비선형성은 더욱 강해지므로 이를 구현하기 위한 완전 비선형 수치수조 개발이 필요하다.

완전 비선형 수조에 대한 연구는 Lunguet-Higgins and Cokelet (1976)의 Mixed Eulerian-Lagrangian 연구 이후, 포텐셜 이론에 기반한 수치연구들이 다양하게 시도되었다. Cointe(1990)는 저차경계요소법을 이용하여 자유수면 유동을 모사하고 채널의 수조시험과 비교하였다. Cao et al.(1991)은 기존 고차요소법과 저차요소법의 요소의 경계적분에 소요되는 시간을 줄이기 위해 특이점을 비특이화시켜 요소적분을 생략하는 방법으로 계산 비용을 줄이는 비특이화 경계적분방정식 방법을 제안하였고 Beck (1994)은 이를 이용하여 완전 비선형 자유수면 유동을 모사하였다. Boo and Kim(1994)은 고차경계요소법을 이용하여 Stokes 3차 불규칙 파에 대한 해석을 수행하였다. Koo and Kim(2004)은 자유수면 또는 고정된 구조물에 대한 해석을 수행한 앞선 제한된 연구들과 달리 자유수면과 부유체의 상호작용을 2차원 완전 비선형 수치 수조에서 저차경계요소법을 이용하여 구현하였다. 국내에서는 Sung et al.(1997)은 고차경계요소법을 이용하여 3차원 슬로싱 문제에 대한 완전 비선형 해석을 수행하고 GMRES (Genralized minimal residual) 알고리즘을 이용하여 선형연립방정식의 효율적인 풀이법을 정립하였다. 그리고 Sung and Choi (2000)은 고차경계요소법을 이용하여 2차원 비선형 방사문제를 모사하고 가속도포텐셜기법을 적용하여 보다 정확하고 안정적인 비선형 방사력을 산정하였다. 이후 Koo and Choi(2009)는 저차경계요소법을 이용하여 2차원 수조의 비선형파의 파형변화와 속도 분포를 해석하여 고찰하였다.

앞선 연구와 같이 경계요소법을 기반한 수치수조는 효율적으로 완전 비선형 수조의 해법으로 널리 사용되고 있다. 경계요소법을 이용할 경우, 완전 비선형 자유수면을 구현하기 위하여 매 시간 간격마다 자유수면을 업데이트하기 때문에 매 시간 간격마다 정의된 경계면에 대한 영향계수행렬을 계산해야하며 이는 많은 계산시간을 필요하게 한다. 하지만 Cao et al.(1991)에 의해 제안된 비특이화 간접경계적분방정식방법(Desingularized indirect boundary integral equation method, DIBIEM)은 기존 경계요소법과 비교하여 특이점을 유체영역 외부에 배치하므로 특이적분이 생략될 뿐만 아니라 모든 면 적분 또는 선 적분 요소가 고립된 점 요소로 치환될 수 있기 때문에 적분계산이 계산절차에서 생략되어 계산비용이 절약된다. 그리고 경계면에서 계산되는 물리량이 다른 이산화된 경계면의 물리량과 연속적이다(Zhang et al., 2006). 특히 고립된 점 요소로 정식화할 경우, 기본해의 형태가 방사기저함수로 정의되기 때문에 무요소 방법으로 치환된다. 무요소 방법의 경우, 요소의 Connectivity에 대한 정보가 필요하지 않기 때문에 경계면 생성과 같은 별도의 연산처리 또한 필요하지 않다. 하지만 완전 비선형 자유표면 경계조건 처리 또는 평활화 기법에 요소의 Connectivity가 필연적으로 필요하다. 비선형 자유표면 경계조건 처리를 위한 자유수면의 구배계산(Wu and Tsay, 2009)과 평활화를 위해 방사기저함수 기법을 도입하면 완전한 무요소 기법의 수치수조를 구현 할 수 있다. 무요소법은 격자생성을 위한 수치기법이 요구되지 않기 때문에 경계면 처리에 더욱 유연한 장점이 있다.

본 연구에서는 비특이화 간접적분방정식방법을 이용하여 2차원 완전 비선형 수치수조를 개발하였다. 앞서 언급한바와 같이 비특이화 간접적분방정식방법은 경계적분 계산의 효율이 높고 수학적으로 단순하기 때문에 코딩을 위한 특별한 기술이 요구되지 않으며 벡터화가 용이하다(Cao and Beck., 2016). 추가로 방사기저함수를 도입하여 격자 재성성이 필요없는 무요소 기반의 수치수조를 구현하였다. 자유수면의 비선형적인 변화를 반영하기 위하여 Mixed Eulerian-Lagrangian(MEL) 정식화를 이용하였다. 시간 적분을 위해 수치 안정성이 뛰어난 Runge-Kutta 4차 방법을 사용하였다. 장시간의 수치계산을 위해 수치 감쇠영역을 적용하여 반사파의 영향을 최소화 하였다. 수치수조의 정확성 검증을 위해 해석 해 및 기 발표된 수리시험에 대한 시뮬레이션을 수행하고 결과들을 비교, 분석하였다.

2. 수학적 정식화 및 수치방법

2.1 경계치 문제

본 연구에서는 Fig. 1과 같이 2차원 직교좌표계(x,z)를 사용하였고 z는 평균수면을 기준으로 상방으로 정의하였다. 자유수면경계, 조파경계, 바닥경계 및 벽 경계조건으로 구성된 수치수조의 경계치 문제를 Fig. 1에서 확인할 수 있다.

Fig. 1.

Schematic diagram for two-dimensional numerical wave tank

유체의 비압축과 비점성 그리고 유동의 비회전을 가정하면 유체속도는 속도 포텐셜의 구배로 표현될 수 있고 지배방정식인 연속 방정식은 식 (1)과 같이 라플라스방정식으로 치환된다.

(1) 2Φ=0influidregion

자유수면 경계면의 압력은 대기압과 동일해야하며 경계면의 속도는 유체입자의 속도와 동일하다. 이러한 물리적 특성을 동역학적 경계조건과 운동학적 경계조건이라 하며 식 (2)(3)과 같다.

(2) Φt=gη12|Φ|2Paρonfreesurfaceboundary
(3) ηt=Φη+Φzonfreesurfaceboundary

여기서 g는 중력가속도, η는 자유표면의 진폭, Pa는 대기압 그리고 ρ는 유체의 밀도를 나타낸다.

수저면과 벽면은 불투과 경계조건으로 식 (4)로 표현된다.

(4) Φn=0onwall&bottomboundary

조파기 경계면은 식 (5)와 같이 조파면 속도 Un으로 정의된다. Piston type 또는 Flap type과 같은 조파기의 형태에 따라 Un을 정의할 수 있다.

(5) Φn=Unonwavemakerboundary

2.2 비특이화 간접경계적분방정식 방법

주어진 경계치 문제를 풀기위해 기본해인 식 (6)과 간접경계 적분방정식 (7)을 이용한다.

(6) G(x,ξ)=12πln(r)wherer=(xξ1)2+(zξ2)2
(7) ϕ(x)=Γσ(ξ)G(x,ξ)dS(ξ)

여기서 소오스 σ(ξ)는 물리적인 의미가 없고 단지 경계조건에 의해 결정되는 미지값이다. 경계면 Γ에서 선 소오스를 비특이화 거리(Disingularized distance) ld로 정의된 특성길이 만큼 유체영역 외부에 분포시키면 기본해인 식 (6)의 특성으로 식 (7)의 선 소오스 적분은 단일 점 소오스의 One node gauss quadrature로 근사 가능하며 식 (8)이 유도 된다.

(8) ϕ(x)=Γσ(ξ)G(x,ξ)dS(ξ)

여기서 σ(ξ)dS의 곱을 σ'(ξ)로 치환하면 식 (8)은 더욱 단순화된 식 (9)가 된다.

(9) ϕ(x)=Γσ(ξ)G(x,ξ)

식 (8) 또는 식 (9)와 같이 특이점이 유체영역 외부에 배치된 경계적분방법을 비특이화 간접경계적분방정식 방법이라 한다. 그리고 식 (8)과 달리 식 (9)의 경우, 기본해가 고립된 점 소오스로 방사기저함수 형태이기 때문에 무요소 방법으로 정의되기도 한다(Wu and Tsay, 2009). 이는 식 (9)의 해법의 경우, 격자가 필요 없으며 경계조건과 소오스의 위치 정보만 요구하기 때문이다. 비특이화된 무요소 계산 경계면은 Fig. 1에서 확인할 수 있다. 추가로 식 (9)의 타당성을 확인하기 위하여 3장에서 수치계산을 수행하고 비교하였다.

식 (9)식 (7)과 비교하여 소오스의 선적분 및 특이점의 계산이 요구되지 않기 때문에 계산 속도가 향상되며 수치코드 작성이 단순해지는 장점이 있다(Cao et al., 1991; Cao and Beck, 2016). 특히 Cao et al.(1991)은 다양한 비특이화 거리(0.8~3.0 (Amesh)1/2)에서 체계적인 연구를 수행하여 수치해석이 수행된 비특이화 거리의 범위에서 평면 패널의 해석적 방법 대비 70~60% 이상의 계산시간이 감소됨을 보고하였다. 그리고 다른 경계면의 근방의 속도가 연속적이고 부드럽게 계산된다(Zhang et al., 2006). 본 연구에서는 ld식 (10)과 같이 정의하였다. 여기서 Lmesh는 국부격자 크기(경계점의 평균거리)를 나타내며 ab는 2차원이므로 모두 1을 이용하였다.

(10) ld=a(Lmesh)b

매 시간 간격마다 식 (9)를 경계조건 식 (2)-(5)에 대입하여 연립방정식을 구성할 수 있다. 구성된 연립방정식을 이용하여 소오스 σ'(ξ)를 결정할 수 있으며 결정된 소오스 σ'(ξ)식 (9)에 대입하여 포텐셜과 경계면에서의 속도를 계산할 수 있다.

2.3 자유수면의 시간적분 방법

완전 비선형 수치수조의 시뮬레이션은 매 시간 간격 마다 자유수면상의 포텐셜과 경계면의 변화를 계산하여야 한다. 본 연구에서는 Mixed Euler-Lagrangian(MEL) 정식화(Lunguet-Higgins and Cokelet, 1976)를 사용하여 자유수면의 경계를 추적하였다. 본 정식화의 경우, Euler 좌표계에서 자유수면 경계조건을 매시간 간격 마다 계산하고 변경된 자유수면의 경계면과 물리량을 Lagrangian 좌표계에서 수정하는 방법이다. 만약 자유수면의 경계면의 속도를 V라 정의할 때 전미분식은 식 (11)과 같이 정의할 수 있다.

(11) δδt=t+V

식 (11)을 이용하여 자유수면 경계조건 식 (2)식 (3)을 변환하면 아래 식 (12)식 (13)이 된다.

(12) δΦδt=gη12ΦΦPaρ+VΦυ(x)Φ
(13) δηδt=(ΦV)η+Φzυ(x)η

여기서 ∇Φ는 입자의 속도를 나타내며 자유수면의 대기압 Pa는 0으로 가정하였다. 본 연구에서는 Fig. 2과 같이 자유수면의 경계면이 수직방향으로만 움직이는 Semi-Lagrangian 정식화를 이용하였고 V식 (14)와 같이 정의할 수 있다. 본 정식화의 경우, 경계면의 노드만 수정하고 요소 Connectivity의 재생성이 요구되지 않는 장점이 있다.

Fig. 2.

Schematic diagram for semi-Lagrangian approach

(14) V=(0,δηδt)

정상상태의 수치시뮬레이션을 지속적으로 수행하기 위해 Fig. 3과 같이 감쇠영역을 설정하였다.

Fig. 3.

Schematic diagram for damping zone

식 (12)(13)에 포함된 v(x)는 감쇠 계수로 식 (15)와 같이 정의할 수 있다.

(15) υ(x)={0,xx0γω(xx0λdamping)2,x>x0

여기서 γ는 조정계수이며 λdamping은 감쇠영역의 길이를 나타낸다. 본 연구에서는 γ는 1을 사용하였고 λdamping은 생성할 파장보다 길게 설정하였다. 시간 적분의 경우, 비선형 자유수면 문제에서 수치안정성이 뛰어난 Runge-Kutta 4차 방법을 적용하였다.

2.4 수치 안정화 기법

안정적인 수치 시뮬레이션을 위해, 조파기의 경계면의 속도를 점진적으로 증가시켜 일정 값에 이르게 하는 램프함수를 식 (16)과 같이 이용하였다.

(16) r(t)={1tTramp12(1cos(πtTramp)t<Tramp

여기서 Tramp은 램프주기를 나타내며 생성하는 파주기보다 긴 주기를 사용하였다. 자유수면의 시뮬레이션의 경우, 자유수면의 경계면의 수치불안정성이 발생되며 이를 억제하기 위한 평활화 방법이 요구된다. 본 연구에서는 B-spline 방사기저함수(Saranli and Baykal, 1998)를 이용한 평활화를 수행하였고 식 (17)와 같이 정리 된다.

(17) ηsmooth(x)=i=1NK(x,xi)η(xi)i=1NK(x,xi)

where

K(x,xi)=14s2{s3+3s2(s|xxi|)+3s(s|xxi|)23(s|xxi|)3,|xxi|s(2s|xxi|)3,s<|xxi|2s0,2s<|xxi|

본 방법의 타당성을 확인하기 위하여 제 1종 베셀함수에 균등 잡음(Uniform noise)을 더한 시험함수에 평활화를 수행하여 Fig. 4에 나타내었다. 시험함수에 대하여 일정 수준의 평활화가 이루어짐을 확인할 수 있다. 본 평탄화 기법을 수치수조에 적용한 결과, 파의 4차항까지 보존됨을 3.3장의 Chapalain et al.(1992)의 수리시험의 계산을 통해 확인하였다.

Fig. 4.

Smoothing filter using B-spline RBF

2.5 자유수면의 구배 계산

Semi-Lagrangian 정식화를 사용할 경우, 자유수면 구배 ∇η의 추가 계산이 필요하다. 기존 방법의 경우 자유수면의 구배를 계산하기 위해 요소의 Connectivity 정보를 이용한 수치 차분이 이용되고 있으며 이는 수치적 부정확성을 증가시키는 것으로 알려져 있다(Beck, 1994). 본 연구에선 이를 극복하기 위하여 방사기저함수(Radial basis function; RBF) 형태의 인공신경망 내삽기법(Wu and Tsay, 2009)을 이용하여 자유수면의 구배를 계산하였다. 자유수면을 내삽할 방사기저함수는 가우시안분포함수를 이용하였다. 본 연구에서 2차원 자유수면의 내삽을 식 (18)과 같이 정의하였다.

(18) η(x)=i=1NFSαifi(x),fi(x)=exp((xxi)2σi2)

여기서 fi(x)x에서 αi의 가중치를 가진 1 차원 방사기저함수이다. Collocation 방법을 통하여 매 시간 간격마다 αi가 계산된다. αi식 (19)을 이용하여 자유수면의 구배를 계산할 수 있다.

(19) ηx(x)=i=1NFSαifi(x)x

본 방법의 타당성을 확인하기 위하여 제 1종 베셀함수의 해석 구배와 수치 구배를 비교하여 Fig. 5에 나타내었다. 식 (18)에서 표준편차 σi를 통해 방사기저함수의 형태와 범위를 조정할 수 있다. 본 연구에서는 σi를 경계점 간격의 2배로 하였다.

Fig. 5.

Numerical derivative using RBF collocation method

3. 수치 계산 및 결과

본 연구에서 개발된 수치 수조의 특성 및 타당성을 확인하기 위하여 수치시험을 수행하고 분석하였다. 우선 비특이화된 선 소오스와 점 소오스 사용의 따른 수치결과를 비교하여 식 (9)의 타당성을 검토하였다. 그리고 Stokes 2차 파에 대한 수치계산을 수행하여 파 진폭을 해석해와 비교하였다. 추가로 조파면의 속도를 Stokes 2차 파의 해석해와 Piston type을 비교하여 공간적 파형 변화의 유무를 확인하였다. 마지막으로 시험 결과가 발표된 Chapalain et al.(1992)의 수리시험 및 Luth et al.(1994)의 수리시험에 대하여 계산을 수행하고 결과를 비교하였다. 관련 수리시험에서 생성된 천해파는 자유수면의 비선형성이 강하기 때문에 수치수조의 완전 비선형 자유수면 특성 검증에 적합하다고 판단된다.

3.1 비특이화된 소오스 종류에 따른 수치계산 비교 및 수치수조의 수렴도 해석

본 연구에서 식 (9)의 타당성을 확인 및 개발된 수치수조의 수렴도 해석을 위해 수치시험을 수행하였다. 수행한 수치 조파 시험 조건은 Table 1과 같으며 h는 수심, Stroke는 조파기 진폭의 2배 그리고 T는 조파기의 주기를 나타낸다. 조파기의 형태는 피스톤 형태로 Fig. 6과 같다.

A simulation condition for a convergence tests

Fig. 6.

The sketch of the numerical wave tank for simulations

우선 2장에서 언급하였던 식 (9)의 타당성을 확인하기 위해 수치시험을 통한 비교를 수행하였다. 식 (8)을 비특이화된 선 소오스(Desingularized line source), 식 (9)를 비특이화된 점 소오스(Desingularized point source)로 정의하고 수치계산을 수행하였다. 조파기로부터 4m 및 8m 떨어진 곳의 파 시계열을 Fig. 7에 도시하였고 두 결과가 일치함으로 식 (9)의 타당성을 확인하였다. 이는 기본해의 특성에서 유추될 수 있는 결과이다. 따라서 본 연구에서 개발된 수치수조는 비특이화된 점 소오스의 커널을 이용하였다. 이는 영향행렬 계산에서 모든 요소의 적분을 생략하고 거리의 함수인 log(r)만 계산하기 때문에 직관적으로 계산량이 줄어듬을 알 수 있다.

Fig. 7.

Time history of wave at specific locations

비특이화거리 ld의 변화에 대한 민감도 해석 역시 수행하였다. 비특이화 거리는 Cao et al.(1991)에서 수행한 범위인 0.8, 1.0 그리고 1.2에 대하여 수행하였고 계산결과는 Fig. 8에 도시하였다. 계산 시 시간스텝 dtT/217인 0.0161s., 파장대비 노드 간격을 λ/61인 0.1m로 고정하였다. 수행된 비특이화 거리 범위 내에서 자유수면은 동일하게 계산되었다. 이는 일정 범위의ld의 변화에 대하여 둔감하다기 발표된 연구(Cao et al., 1991)와 동일하다.

Fig. 8.

Sensitivity test of desingularized distance

주요 변수의 수렴도 해석을 위해 파장대비 노드수는 61, 51, 34 및 17에 대하여 변화를 주어 계산을 수행하였고 그 결과를 Fig. 9에 도시하였다. 계산 시 시간스텝 dtT/217인 0.0161s로 고정하였다. 파장대비 노드수가 증가함에 따라, 61개 노드를 기준으로 17개 노드에서 8.24% 차이가 그리고 51개 노드에서 0.68% 차이로 수렴하고 있음을 확인 할 수 있다.

Fig. 9.

Covergence study according to the variation of number of node per wave length

주기대비 시간스템에 대하여 217, 87 및 43에 대하여 변화를 주어 계산을 수행하였고 그 결과를 Fig. 10에 도시하였다. 계산 시 파장대비 노드 간격을 λ/61인 0.1m로 고정하였다. 주기대비 시간간격이 작아질수록 수렴함을 확인할 수 있다.

Fig. 10.

Covergence study according to the variation of time step per period

3.2 Stokes 2차 파의 비교 수치계산

개발된 2차원 수치수조의 검증을 위해 해석해가 존재하는 Stokes 2차 파에 대한 수치계산을 수행하였다. 계산에 사용된 조건은 Goda(1998)에 의해 수행된 수리시험 조건을 사용하였고 Table 2에 정리하였다. 조파면의 속도는 Stokes 2차 파 이론을 이용하여 식 (19)와 같이 설정하였다. 조파면에서 4m 떨어진 곳의 계산된 파 진폭 시계열을 선형 파 및 Stokes 2차 파와 함께 Fig. 11에 도시하였다.

A simulation condition for 2nd order Stokes wave

Fig. 11.

Time history of wave at specific location (4 m from wave maker)

(19) Un=Φx=gkH2ωcosh(k(h+z))cosh(kh)cos(kxωt)+3ωkH216cosh(2k(h+z))sinh4(kh)cos(2(kxωt))

수치수조의 생성된 파의 주기는 해석해와 전반적으로 잘 일치하고 파의 형상은 특히 Stokes 3차 파와 비교하여 파정과 파저가 아주 유사함을 확인할 수 있다.

일반적으로 조파수조의 조파면의 속도는 Stokes 2차 파의 이론을 구현하기 힘든 Piston type 또는 Flap type이 대부분이다. Goda(1998)Koo and Choi(2009)의 의하면 시험 수조에서 생성되는 비선형 천해파의 경우, 조파기 주변의 입자의 비선형 변위와 조파기면의 변위 불일치로 인해 2차항의 자유파가 발생하기 때문에 수조 길이방향으로 공간적 파형변화가 발생한다. Fig. 12는 조파면의 속도를 Stokes 2차 파의 이론해와 Piston type의 속도를 이용하여 시뮬레이션하고 수조 길이방향으로 조화진폭을 계산하였다. Piston type의 조파면 속도를 적용한 경우, 공간적 파형 변화가 1~2차항 모두 발생하며 Goda(1998)의 시험과와 Koo and Choi(2009)의 계산과 잘 일치함을 확인하였다. Stokes 2차 파의 이론해를 조파면 속도에 적용한 경우, 공간적 파형 변화가 사라짐을 확인하였고 이는 Goda(1998)Koo and Choi (2009)의 해석과 일치한다.

Fig. 12.

Comparison of spatial variation of harmonic amplitude

3.3 Chapalain et al.(1992)의 수리시험

Chapalain et al.(1992) 수리 시험을 Fig. 13과 같이 수치수조로 구현하였다. 조파수조의 천해파의 시험의 경우, 비선형성이 강하기 때문에 공간적 파형변화가 두드러진다. 따라서 Chapalain et al.(1992) 수리 시험은 완전 비선형 수치수조의 비선형 재현성을 검증하기 적합하다. 본 연구에서 수행된 수치계산의 시험 조건은 Table 3과 같다. 여기서 h, eT는 수심, 조파기 변위의 진폭 그리고 조파기 변위의 주기를 나타낸다. 본 시뮬레이션의 경우, 피스톤 형태의 조파기를 구현하였다. Fig. 14은 공간적 파형변화에 대한 계산결과와 수리 시험결과를 함께 도시하였으며 4차 조화항까지 잘 일치함을 확인하였다.

Fig. 13.

The sketch of the numerical wave tank for Chapalain’s experiment(1992)

Chapalain’s experimental characteristics(1992)

Fig. 14.

Spatial variation of harmonic amplitude

3.4 Luth et al.(1994)의 수리시험

Luth et al.(1994)의 수리 시험을 수치수조로 구현하여 Fig. 15와 같이 나타내었다. Submerged bar의 제원은 높이 0.3m 그리고 진입경사면이 1:20의 비율과 유출경사면이 1:10의 비율을 가진다. Luth et al.(1994)의 수리시험과 동일한 파 진폭 0.01m와 파 주기 2.02sec.의 조건에서 수치 계산을 수행하였다. 수리시험과 동일하게 Piston type의 조파기로 시뮬레이션을 수행하였고 조파기 이론을 통해 조파기의 Stroke는 0.29859m를 결정하였다. 수리시험과 비교하기 위하여 조파면에서 4m, 12.5m, 14.5m, 17.3m과 21.0m 떨어진 지점의 파진폭 시계열을 비교하여 Fig. 16에 도시 하였다. 조파면에서 4m 떨어진 지점은 Reference point로 수면 변화가 없기 때문에 조파기 이론과 동일한 0.01m의 비선형 파를 생성함을 확인하였고 수리시험과 유사한 결과를 나타냄을 확인하였다. 다른 지점의 파진폭 시계열 또한 Luth et al.(1994)의 수리시험과 잘 일치함을 확인할 수 있다. 특히 조파면에서 12.5m와 14.5m 떨어진 지점에서는 수심이 낮아져서 파고가 증가하고 파장이 짧아지며 17.3m과 21.0m 떨어진 지점에서는 수심이 깊어져 파고가 낮아지는 현상을 수치시계산에서도 관찰할 수 있었다.

Fig. 15.

The sketch of the numerical wave tank for Luth’s experiment (1994)

Fig. 16.

Wave elevations at specific locations

4. 결 론

본 논문에서는 DIBIEM을 이용하여 2차원 완전 비선형 수치수조를 개발하였다. 자유수면의 구배계산과 안정성을 위한 평활화 방법에 방사기저 함수를 사용하였다. 개발된 수치기법은 무요소 기반의 방법으로 분류되며 계산을 위한 격자가 필요 없기 때문에 기존 방법과 비교하여 수치적으로 유연하다.

수치수조의 특성 파악 및 검증을 위해 다양한 수치계산들을 수행하였다. 우선 비특이화된 소오스의 종류에 따른 수치결과를 비교하였고 수치적으로 계산이 단순한 점 소오스가 선 소오스를 대체할 수 있음을 확인하였다. 그리고 Stokes 2차 파 이론을 이용하여 수치수조의 시계열을 해석해와 비교 검증하였고 조파면의 속도를 Stokes 2차 파 이론해를 이용하면 Piston type의 조파기를 사용할 때 발생하는 공간적 파형변화가 일어나지 않음을 확인하였다. 마지막으로 비선형성이 강한 천해역의 수리시험인 Chapalain et al.(1992)의 수리 시험 및 Luth et al.(1994)의 수리시험에 대하여 수치계산을 수행하였다. 계산결과들을 기 발표된 수리시험결과와 비교하였으며 정량적으로 잘 일치함을 확인하였다.

추후 본 연구에서 성공적으로 개발된 방사기저 함수 기반의 2차원 완전 비선형 수치수조를 이용하여 심해공학수조시험의 설계 및 해석에 활용하고 3차원 수치수조로 확장연구를 수행할 예정이다.

Acknowledgements

본 연구는 해양수산부의 국가개발사업인 ‘심해공학수조수조 운영을 위한 연구 인프라 구축 및 심해플랜트 Pre-FEED 원천핵심기술개발’ (PMS3850)의 결과물임을 밝히는 바입니다.

References

Beck RF. 1994;Time-domain Computations for Floating Bodies. Applied Ocean Research 16:267–282.
Boo SY, Kim CH. 1994. Nonlinear Irregular Waves and Forces on Truncated Vertical Cylinder in a Numerical Wave Tank. In : Proceedings of 7th International Offshore and Polar Engineering Conference. Honolulu, HI. ISOPE; 3p. 76–84.
Cao Y, Beck RF. 2016;Desingularized Boundary Integral Equations and Their Applications in Wave Dynamics and Wave-body Interaction Problems. Journal of Ocean Engineering and Science 1(1):11–29.
Cao Y, Schultz WW, Beck RF. 1991;Three Dimensional Desingularized Doundary Integral Methods for Potential Problems. International Journal for Numerical Methods in Fluids 12:785–803.
Chapalain G, Cointe R, Temperville A. 1992;Observed and Modeled Resonantly Interacting Progressive Water-waves. Coastal Engineering 16:267–300.
Cointe R. 1990;Numerical Simulation of a Wave Channel. Engineering Analysis with Boundary Elements 7(4):167–177.
Ducrozet G, Bonnefoy F, Touze DL, Ferrant P. 2012;A Modified High-Order Spectral Method for Wavemaker Modeling in a Numerical Wave Tank. European Journal of Mechanics B/Fluids 34:19–34.
Goda Y. 1998. Perturbation Analysis of Nonlinear Wave Interactions in Relatively Shallow Water. In : Third International Conference on Hydrodynamics. Seoul, South Korea. p. 33–51.
Koo WC, Kim MH. 2004;Freely Floating-body Simulation by a 2D Fully Nonlinear Numerical Wave Tank. Ocean Engineering 31:2011–2046.
Koo WC, Choi KR. 2009;Spatial Modulation of Nonlinear Waves and Their Kinematics using a Numerical Wave Tank. Journal of Ocean Engineering and Technology 23(6):12–16.
Longuet-Higgins MS, Cokelet ED. 1976;The Deformation of Steep Surface Waves on Water: I. a Numerical Method of Computation. Proceedings of the Royal Society of London, Series A 350:1–26.
Luth HR, Klopman G, Kitou N. 1994. Kinematics of Waves Breaking Partially on an Offshore Bar; LDV Measurements of Waves with and without a Net Onshore Current. Report H-1573. Delft Hydraulics;
Wu NJ, Tsay TK. 2009;Applicability of the Method of Fundamental Solutions to 3-D Wave-body Interaction with Fully Nonlinear Free Surface. Journal of Engineering Mathematics 63:61–78.
Saranli A, Baykal B. 1998;Complexity Reduction in Radial Basis Function(RBF) Networks by Using Radial B-spline Functions. Neurocomputing 18:183–194.
Sung HG, Hong SY, Choi HS. 1997;A Study on the Boundary Element Method for Numerical Analysis of Nonlinear Free Surface Waves(I). Journal of the Society of Naval Architects of Korea 34(4):53–60.
Sung HG, Choi HS. 2000;Numerical Analysis of Two-Dimensional Nonlinear Radiation Problem Using Higher-Order Boundary Element Method. Journal of the Society of Naval Architects of Korea 37(1):67–81.
Zhang XT, Khoo BC, Lou J. 2006;Wave Propagation in a Fully Nonlinear Numerical Wave Tank: A Desingularized Method. Ocean Engineering 33:2310–2331.

Article information Continued

Fig. 1.

Schematic diagram for two-dimensional numerical wave tank

Fig. 2.

Schematic diagram for semi-Lagrangian approach

Fig. 3.

Schematic diagram for damping zone

Fig. 4.

Smoothing filter using B-spline RBF

Fig. 5.

Numerical derivative using RBF collocation method

Fig. 6.

The sketch of the numerical wave tank for simulations

Fig. 7.

Time history of wave at specific locations

Fig. 8.

Sensitivity test of desingularized distance

Fig. 9.

Covergence study according to the variation of number of node per wave length

Fig. 10.

Covergence study according to the variation of time step per period

Fig. 11.

Time history of wave at specific location (4 m from wave maker)

Fig. 12.

Comparison of spatial variation of harmonic amplitude

Fig. 13.

The sketch of the numerical wave tank for Chapalain’s experiment(1992)

Fig. 14.

Spatial variation of harmonic amplitude

Fig. 15.

The sketch of the numerical wave tank for Luth’s experiment (1994)

Fig. 16.

Wave elevations at specific locations

Table 1.

A simulation condition for a convergence tests

Item Value
h [m] 0.4
Stroke [m] 0.226
T [s] 3.5
Wave maker type Piston type

Table 2.

A simulation condition for 2nd order Stokes wave

Item Goda’s Experiment (1998)
h [m] 0.25
Amplitude [m] 0.025
T [s] 3.5

Table 3.

Chapalain’s experimental characteristics(1992)

Trial h [m] e [m] T [s]
A 0.4 0.078 2.5
C 0.4 0.113 3.5
D 0.3 0.078 2.5
H 0.4 0.078 3.0