2.1 파랑모형

ihFoam 파랑모형인 RANS, 질량 보존식, VOF식을 기술하면 다음과 같다.

여기서 U는 속도 벡터, g는 중력가속도, ρ는 밀도, μ_eff = μ + ρν_turb는 동점성계수, p^*는 유사 동적압력, X는 위치벡터, α₁는 각각의 단위 Cell에서 해수가 점유하고 있는 부분을 지칭하는 VOF 함수를 각각 나타낸다.

2.2 Porous media에서의 파랑모형[VARANS equations in iH-3VOF]

항 외곽시설과 잠제, 돌제, 이안제 등과 같은 연안 보호공의 상당부분은 사석층으로 구성되며, 이 경우 다공성 매질에서 발생하는 형상항력과 관성항력으로 인해 파랑모형은 수정되어야한다. Darcy(1856)의 사질토 내부에서의 유체 흐름에 관한 연구가 알려진 후, 다공성 매질에서의 유체 해석을 형상항력으로 기술하려는 흐름은 이후 다공성 매질에서의 유동해석에 큰 흐름을 형성한다. 층류의 경우 선형 형상항력만으로도 해석이 가능하나, Forcheimer(1901)는 큰 Reynolds수에서도 적용이 가능하도록 이차 형상 항력항을 추가하였으며, Polubarinova와 Kochina(1962)는 1901년에 Forcheimer가 제시한 식을 비정상류로 확장하기 위해 관성력을 추가하였다.

본고에서는 del Jesus et al.(2012)에 의해 제시된 VARANS[Volume Averaged Reynolds Averaged Navier Stokes Eq.]를 사용하였으며 이를 기술하면 다음과 같다.

여기서 u는 평균유속 또는 Darcy 속도, n은 공극율, A와 B는 형상항력 계수, α₁은 VOF 계수, C는 관성력 형상항력계수로 물성에 종속한다.

다공성 물질에서 비어있는 공간에 존재하는 유체의 속도를 의미하는 실제속도 U는 평균유속 또는 Darcy 속도를 나타내는 u와 다음과 같은 관계식을 추종하며

선형 형상항력계수 A, 비선형 형상항력 계수 B는 다음과 같이 정의되며

여기서 D₅₀은 다공성 물질의 명목 직경, KC는 Keulegan-Carpenter수를 각각 나타내며 KC는 다음과 같이 정의되며

, u_M은 최대 진동 속도, T_o는 진동주기를 나타내며, KC는 vortex induced vibration으로 인한 추가적인 마찰을 고려하기 위해 도입되었다.

마찰 매개변수들과 연관되어 α와 β는 조정되며, 관성력 계수 C값은 0.34로 유지된다.

Fig. 2에는 이해를 돕기 위해 d₅₀에 따른 항력계수 A, B의 변화를 표시하였다.

Fig. 2.

Trajectories of damping coefficient A and B by varying d₅₀ for specific void ratio.

본고에서는 del Jesus et al.(2011)의 실험결과를 토대로, 방파제를 구성하는 사석의 명목직경 D_n50 = (M₅₀/ρ_r)^1/3은 0.015 m, 투수계수는 0.41, 선형 형상항력계수 A = 10,000, 비선형 형상항력 계수 B = 6로 취하였으며, 여기서 M₅₀는 누적분포에서 중앙에 해당되는 사석 질량, ρ_r는 사석 밀도를 각각 나타낸다.

2.3 Energy absorbing Boundary condition

계산영역내부에서 반사되어 외해로 회귀하는 반사파가 개방경계에서 적절히 제어되지 않을 경우 수치모의 결과는 왜곡될 수 있다.

반사파를 제어하기 위한 방법으로는 passive wave absorbers와 active wave absorption 등으로 구분된다.

Artificial beach, 다공성 물질, 다공성 판, 수치적으로 감소시키는 방법 등이 passive wave absorbers에 해당한다. 이러한 방법은 완벽히 파랑 에너지를 소멸시키지 못하며 반사파가 발생할 수 있다. 특히, 장주기파의 경우 artificial beach로 는 완벽히 소산되지 못해 반사파가 생성되는 것으로 알려진바 있다(Pengzhi and Liu, 1999; Lara et al., 2006; Wei and Kirby, 1995; Losada et al., 2008).

Active wave absorption(Schaffer and Klopman, 2000; Torres et al., 2010; Lara et al., 2011)은 조파 피스톤의 움직임을 측정된 강도에 맞게 조정하여 입사파의 재반사를 제어하는 것으로, 처음에는 수리모형실험을 위해 개발되었으나, 수치모형에서도 조파 경계와 개방 경계에서도 적용가능하며, 경계에 적절한 유속성분을 수치적으로 부과하여, 파랑을 흡수할 수 있다.

Active wave absorption systems 유도과정을 살펴보면 다음과 같다.

먼저 다음과 같이 정의되는 depth averaged mass conservation Eq.과

ζ = Ae^i(kx−ωt), U~ = U_oe^i(kx−ωt)로부터 다음과 같은 관계식을 얻을 수 있으며, 여기서 A는 진폭, U_o는 유속 진폭, k는 파수, ω는 각 주파수를 나타낸다.

상기식에서 U~ 는 depth averaged 유속, h는 수심을 나타내며, h와 ζ는 관측치로 구한 기지 값, U~ 는 미지 값으로 취급한다.

파속 c는 선형이론으로부터 다음과 값이 기술할 수 있으나

여기서 파수 k는 관측치로부터 추정하기가 매우 어려우나, 최근 Wellens(2012)은 상기식을 다음과 같이 근사하여 전술한 어려움을 해결하였다.

상기식에서, c는 천해 지역에서 gh로 수렴한다.

반사파를 억제하기 위해서 경계에서는 입사하는 파랑의 반대방향으로 입사 파랑의 속도와 같게 조파를 시킨다. 식 (12)를 반사 파랑 ζ_R을 중심으로 정리하는 경우 다음과 같으며

, 여기서 U_c는 경계에서 수직하게 계산영역 내부로 향하는 방향으로의 correction velocity를 나타내며, 반사파 파고 ζ_R는 다음과 같이 산출할 수 있다.

여기서 ζ_M는 관측 파고, ζ_Y는 목표파고를 나타낸다.

그러나 대부분의 경우 반사파랑은 개방경계 혹은 조판경계로 직각으로 접근하는 것은 아니므로, Schaffer와 Klopman(2000)은 상기식을 다음과 같이 일반화 하였으며,

여기서 ∇β는 반사 파랑의 입사각을 의미하나, 대부분의 경우, 반사 파랑의 입사각도를 미리 예측할 수 없다는 단점을 지닌다[Fig. 3 참조].

Fig. 3.

Schematic sketch of the corrected velocity at the wave paddle or open boundary.

이러한 어려움은 각 조파판 전면수역에서 자유수면을 측정하고, 이를 토대로 digital filter를 사용하여 파랑의 진행방향을 산출한다.

이렇게 추정된 방향을 사용하여, 유속을 조파판 직각방향과 접하는 방향으로 다음과 같이 분해하면 파랑의 방향성은 다음과 같이 나타낼 수 있다.

여기서 Ucalc는 총 예상 유속의 절대 값을 나타낸다. 이 값을 두 가지 수평성분으로 나누면 다음과 같다.

여기서 U_corr 값을 제외하고는 기지 값이며, U_corr 값은 파랑을 흡수하기 위해 패달에 수직한 방향으로 적용되는 correction velocity으로 다음과 같이 정의된다.

3.1 수치모형의 검증

2장에서 소개한 수치모형을 검증하기 위해 파향제어식 해안보호공의 주 구조물로 검토되고 있는 직립식 이안제와 잠제를 투과하는 파랑을 수치모의 하였으며, 이와 더불어 수리모형실험도 병행하여 수행하였다. 직립식 이안제와 잠제는 일정 수심부에 거치하였으며[Fig. 4 참조], 수리모형실험에서는 직접적인 조파가 어려운 상대적으로 큰 파고를 유도하기 위해 기울기가 1:30인 사면부를 구조물 전면 수역에 배치하였다. Fig. 5에는 수치모의된 연안보호공 배후수역[Gauge 6, Fig. 4 참조]에서의 해수면 변위를 도시하였으며, 비교를 위해 수리모형실험에서 관측된 시계열자료도 함께 도시하였다.

Fig. 4.

Schematic sketch of the numerical wave flumes.

Fig. 5.

Comparison of the numerically simulated transmitted waves with the measured one in the hydraulic model test [green solid line: hydraulic model test, blue solid line: openFoam].

파고가 상대적으로 적은 경우 연안보호공의 배후 수역에서 관측되는 투과파는 상당히 정확히 재현되는 것을 확인할 수 있다[Fig. 5(a), (b), (d), (e)]. 파고가 상대적으로 큰 경우에도 비교적 정확히 재현되었으나, 미세한 위상차와 파고차가 발생하는 것을 알 수 있다. 이와 더불어 광폭 잠제가 거치된 경우 수리모형 결과에는 고차 조화성분이 관측되나, 수치모의 결과에서 고차 조화성분은 찾아볼 수 없다[Fig. 5(g), (h)].

이렇게 상이한 파형과 파고는 비선형 천수과정에서 생성된 bound mode의 고차 조화성분이 쇄파로 인해 구속이 해제되어 free mode로 전환되는 경우 발생하는 것으로 보이며 이를 정리하면 다음과 같다.

수치모의는 선형이론에 의해 조파된 파랑을 다루는 반면에 수리모형에서는 기울기가 1:30인 일정 사면부가 증폭을 위해 거치되어 있다. 이 경우 사면부에서의 비선형 천수과정과 이 과정에서 진행되는 성분 파랑간의 상호작용으로 bound mode의 고차 조화성분이 파동계에 출현한다. 이렇게 생성된 bound mode의 고차 조화성분은 잠제로 인해 쇄파되는 경우 구속이 해제되어 free mode로 전환되며, 이 경우 파동계를 구성하는 각각의 파랑은 위상차를 지니며, 이는 상쇄 간섭으로 이어져 파고는 감소한다. 이러한 추론은 고차 조화성분이 확연히 관측되는 Fig. 5(g), (h)에서 확인할 수 있다.

전술한 위상차를 설명하기 위해 잠제 혹은 이안제에서의 damping이 반영된 분산관계식을 새로이 유도하였으며 이는 3.2절에 자세히 다루나, 논의를 위해 유도결과를 미리 전하면 파속은 Damping으로 인해 A'sinτ/2kgh만큼 감속되며 감속되는 양은 파형과 마찰력 간의 위상차 τ에 비례한다[Fig. 6 참조].

Fig. 6.

Definition sketch of the phase lag between the damping friction and free surface elevation.

3.2 Porous 층에서의 damping이 파속에 미치는 영향

Porous 층에서의 Damping으로 인한 파속변화 여부를 살펴보기 위해 Damping을 반영한 천수방정식의 분산관계식을 새로이 유도하였으며 이를 정리하면 다음과 같다.

방파제가 천해역에 거치된다는 점을 고려하면 파랑모형은 다음과 같이 간결하게 근사할 수 있다.

여기서 U~ 는 수심평균 유속으로 다음과 같이 정의되며

A와 B는 각각 선형 형상항력계수, 비선형 형상항력 계수를 나타내며, 식 (8), (9)에 기정의 한 바 있다. 해석의 편의를 위해 비선형 형상항력 항을 등가 선형화 기법[equivalent linearization]을 사용하여 선형화하는 경우 식 (22)는 다음과 같이 기술될 수 있으며,

여기서 A' = A + 8B/3π.

식 (21), (24)에서 일방향 진행파라 가정하면 해수면 변위 ζ와 유속 U~ 는 다음과 같이 기술할 수 있다.

여기서 식 (25), (26)을 식 (21), (24)에 대입하고 정리하면 다음과 같은 관계식을 얻을 수 있다.

상기한 계에서 해가 존재하기 위해서는 행렬의 determinant가 영이 되어야 한다는 Solvability condition으로부터 다음과 같은 분산 관계식을 얻을 수 있다.

식 (28)에서 O(1)에서 ω≈kgh이므로 O(ε = A_D/h)에서 ω²는 다음과 같이 근사할 수 있다.

식 (29)의 양변을 ω_o²=k²gh으로 나누어 정규화하고, 섭동법을 사용하여 근사하면 다음과 같은 관계식을 얻을 수 있다.

여기서 파형과 damping 항간의 위상차 e^−iτ[Fig. 6 참조]를 고려하는 경우, 상기식은 다음과 같이 전환된다.

식 (31)에 기초하여 주파수 ω를 실수부와 허수부로 분리하여 다음과 같이 기술하고

여기서

식 (32)을 식 (25)에 대입하고 정리하면 다음과 같은 관계식을 얻을 수 있다.

따라서 식 (32)의 허수부는 Porous층에서의 Damping으로 인한 감쇄를 의미하며, 파속 c는 다음과 같이 수정된다.

여기서 c_o = ω_o/k이며, 따라서 Damping으로 인해 파속은 A'sinτ/2k gh만큼 감속되는 것을 알 수 있다.

3.3 잠제 폭이 투과율과 set-up 량에 미치는 영향

상기한 검증과정을 통해 최적화된 수치모형을 활용하여 파향 제어식 연안보호공의 주 구조물인 잠제 인근 수역에서의 수리특성을 확인하기 위한 수치모의를 수행하였으며, Fig. 7에는 수치모의에서 얻은 지유수면 snapshots을 도시하였다.

Fig. 7.

Snapshots of numerically simulated wave field over the submerged wave breaker.

Fig. 8에는 잠제폭에 따른 투과파 변화를 도시하였으며, 잠제 폭이 증가할수록 감소하는 투과율을 정량적으로 확인할 수 있다.

Fig. 8.

Comparison of transmitted waves by varying the width of submerged breakwater [blue solid line: B = 0.2 m, green solid line: B = 0.6 m, red solid line: B = 1.0 m, green solid line: pale blue-green solid line: B = 1.6 m ].

Fig. 9에는 잠제폭에 따른 잠제 배후수역에서의 수위 상승량 변화[set-up]를 도시하였다. 잠제가 쇄파역 내측에 거치되는 경우 처오름 정점에서 회귀하는 흐름은 이안수로로 집중되어 이안수로 바깥쪽에는 자연해안에서 관측되는 저류는 형성되지 못한다. 이렇게 사라진 저류로 인해 잠제를 가로지르는 방향에서 주기 평균된 유동계의 유속방향은 연안을 향하며[Stokes drift], 종국에는 잠제를 월류하여 배후지로 유입되며, 이는 배후수역 수위 상승으로 이어지는 것으로 판단된다[Fig. 9(a) 참조].

Fig. 9.

Definition sketches of under-tow and variation of set-up at the down wave side of submerged breakwater [blue circle: h_s = 0.14 m, gold circle: h_s =0.17m].

3.4 왜도된 비선형성 파랑이 투과율에 미치는 영향

왜도된 비선형성 파랑이 투과율에 미치는 영향을 해석하기 위한 수치모의를 수행하였다. 현재 transfer function을 이용한 왜도된 파랑의 직접적인 조파는 불가능하며, 따라서 수치모의는 1:30의 기울기를 지니는 천퇴를 포함한 수치수조에서 수행하였다[Fig. 10(a) 참조]. 수치모의는 먼저, 왜도된 비선형성 파랑의 조파여부를 확인하기 위해 방파제가 거치되지 않은 경우에 대해 수행하였다[RUN1]. 이어 왜도된 비선형성 파랑이 투과율에 미치는 영향을 해석하기 위해 동일한 수치수조에서 방파제가 거치된 경우에 대해 수치모의를 수행하였다[RUN2]. 또한 비교를 위해 RUN1에서 방파제 거치예정 수역에서 관측된 파고를 조파파고로 활용하여 수심이 일정한 수치수조를 대상으로 한 수치모의도 병행하여 수행하였다[RUN3, Fig. 10(b) 참조].

Fig. 10.

Schematic sketch of the numerical wave flumes.

RUN1과 2의 경우 입사파고와 주기는 각각 H = 0.11, 0.13, 0.15 m, T = 1.7, 2.2초로 선정하였으며, 파형은 수심[0.57 m]을 고려하여 Stokes I type으로 취하였다. RUN3의 경우에는 보다 정확한 비교를 위해 수심[0.17 m]을 고려하여 천해역에서의 대표적인 비선형 파랑 모형이나 왜도가 영인 Cnoidal wave로 취하였다. Korteweg-de Vries식의 해석해로 정의되는 Cnoidal wave를 기술하면 다음과 같다(Mei, 1989).

여기서 ζ는 해수면 변위, H는 파고, c는 파속, L은 파장, T는 주기, m은 0과 1 사이의 값을 지니는 elliptic parameter, K_m은 제1종 완전 타원적분, E_m은 제2종 타원적분을 각각 나타내며 이를 기술하면 다음과 같으며,

K_m과 E_m은 값에 종속한다.

식 (37)에서 Cn[χ; m]은 Jacobi elliptic function을 나타내며 다음과 같이 정의된다.

여기서 χ와 Ψ는 다음과 같은 관계식을 충족하며

Cnoidal wave 파형은 m이 0에 수렴하면 Stokes I 파형에, m이 1에 수렴하면 Solitary wave 파형에 수렴하며, Cnoidal wave에서의 유속은 다음과 같이 기술될 수 있다.

여기서 ζ_x, ζ_xx, ζ_xxx는 x에 대한 자유수면의 1차, 2차, 3차 derivatives를 각각 나타내며, ∙는 파랑 한 주기 동안의 평균 연산자를 나타낸다.

Fig. 11에는 방파제가 미 거치된 경우[RUN1_A, B, C]에 대해 수행한 세 번의 수치모의에서 얻은 snapshots을 도시하였다. 목표한 파고가 정확히 조파되어 정연하게 이행되는 것을 알 수 있다. 또한, 수치수조 down wave side의 개방경계에 적용된 Energy absorbing Boundary condition가 정확히 구현되어 접근하는 파랑이 수치적 반사 없이 자연스럽게 계산영역을 빠져나가는 것을 확인할 수 있다.

Fig. 11.

Snapshots of numerically simulated wave field.

Fig. 12에는 RUN1_A, B, C, D, E, F에서 방파제 위치[Gauge 14, Fig. 10 참조]에서 관측된 수위 시계열 자료를 도시하였다. 우리의 예상대로 상당히 peaky하고 왜도된 파형이 관측된다. Hoefel과 Elgar (2003)에 준거하면, 정규화된 무차원 파랑 왜곡도[Skew]와 왜도[Assy]는 각각 다음과 같이 기술할 수 있으며

Fig. 12.

Sample time series of the numerically simulated water surface displacement at Gauge 14.

여기서 <•>는 Ensemble mean, ∙~ 는 Hilbert transform을 각각 나타내며, ζ_RMS는 다음과 같이 정의된다.

산출된 skewness와 asymmetry는 Table 2에 정리하였다.

Table 2.

List of the measured skewness and asymmetry of the numerically simulated wave field at Gauge 14

Fig. 13에는 방파제 위치[Gauge 14]에서 관측된 파고와 목표한 파고 사이의 회귀 분석결과를 도시하였다. 1:30의 경사부에서 진행되는 쇄파로 파고가 감쇄되는 것을 알 수 있으며 회귀 분석결과는 다음과 같다.

Fig. 13.

Results of the regression analysis between the measured wave heights and target wave heights [RUN1_A, B, C].

여기서 m는 관측파고, t는 목표 파고를 각각 나타낸다.

Fig. 14에는 t = 40 s에서 RUN1_C, RUN2_C, RUN3_C에서 관측된 파랑 진행방향 유속성분의 등고선도를 도시하였다. 천수과정에서 증가한 강비선형성으로 인해 저면까지 상당한 유속이 유지되는 것을 알 수 있다. 또한 동일한 파고일 경우 Skewed waves의 경우에서 Cnoidal waves의 경우보다 상대적으로 빠른 유속이 관측된다.

Fig. 14.

Contour plots of the horizontal components of numerically simulated wave induced velocities in RUN1_C, RUN2_C, and RUN3_C.

Fig. 15에는 RUN1_A, B, C에서 관측된 평균해면의 변화를 경사부를 중심으로 도시하였다. Radiation stress의 변화에 따라 평균해면이 set-down되고 set-up되는 것을 확인할 수 있으며, 이러한 변화는 파고가 증가할수록 커지는 것을 알 수 있다.

Fig. 15.

Variation of mean water level [MWL] due to set-down and set-up [RUN1_A, B, C, see Table 1].

Table 1.

List of wave conditions used in the numerical simulation

Fig. 16에는 RUN2_A, B, C, D, E, F와 RUN3_A, B, C, D, E, F에서 방파제 배후 수역[Gauge 16]에서 관측된 자유수면 시계열 자료를 도시하였다.

Fig. 16.

Comparison of the transmitted waves at the down-wave side of breaker by skewed waves with the one by the corresponding Cnoidal waves.

Fig. 17, 18, 19, 20에는 파고가 투과율에 미치는 영향을 확인하기 위해 RUN2_A, B, C, RUN3_A, B, C에서 방파제 배후 수역에서 관측된 자유수면 시계열 자료를, Fig. 21에는 투과율[K_t = H_r/H_i]을 각각 도시하였다.

Fig. 17.

Comparison of the transmitted waves at the down-wave side of breaker by skewed waves for varying wave heights [RUN2_A, B, C].

Fig. 18.

Comparison of the transmitted waves at the down-wave side of breaker by Cnoidal waves of varying heights [RUN3_A, B, C].

Fig. 19.

Comparison of the transmitted waves at the down-wave side of breaker by skewed waves of varying heights [RUN2_D, E, F].

Fig. 20.

Comparison of the transmitted waves at the down-wave side of breaker by Cnoidal waves for varying wave heights [RUN3_D, E, F].

Fig. 21.

Comparison of the transmission rate of breaker by Cnoidal waves with the one by Skewed waves for varying wave heights.

파고가 클수록 파속은 증가하는 것으로 보이며, 이는 진폭분산에 기인한다[식 (32), Fig. 17, 18, 19, 20 참조]. Skewed waves가 Cnoidal wave에 비해 파속이 느려지는 현상은 Skewed waves의 경우 Cnoidal wave에 비해 파형과 마찰력간의 위상차가 상대적으로 커서 발생하는 것으로 보인다[3.2절 참조].

그 동안 추정되어왔듯이, T = 1.7 s의 경우 왜도된 비선형파랑이 큰 관성력으로 인해 상대적으로 큰 투과율을 보이는 것을 확인할 수 있으나, T = 2.2 s의 경우 H = 0.15 m를 제외하고는 Cnoidal waves가 Skewed waves에 비해 상대적으로 큰 투과율을 보였다. 이렇게 상이한 결과를 야기한 수리적 배경을 설명하면 다음과 같다.

방파제 전면 수역에서 Cnoidal waves와 Skewed waves가 동일한 파고를 지니도록 조절했으나, 천수과정을 통해 생성된 Skewed waves의 경우 작지 않은 set-up 양[Fig. 15 참조]에서 알 수 있듯 쇄파과정을 거친 것으로 보인다. 이는 Skewed waves의 경우 천수과정에서 파동계에 생성된 bound mode의 고차 조화성분이 free mode로 전환되었음을 의미하나, Cnoidal waves를 구성하는 고차 조화성분은 여전히 bound mode를 유지하며 진행된다. 전술한 현상은 이렇게 생성된 free mode의 고차 조화성분간에 진행되는 상쇄 간섭에 기인하는 것으로 추정되며, 이러한 가설과 주기가 길수록 천수과정이 상대적으로 일찍 시작된다는 사실에 기초하면 T = 1.7 s와는 다르게 T = 2.2 s에서 Cnoidal waves가 상대적으로 큰 투과율을 보이는 현상을 설명할 수 있다.