1. 서 론

일반적인 사석방파제 또는 사석호안은 제체내부가 공극매체로 이루어져 있으며 입사파 에너지가 제체내로 침투하여 내부에서 반사, 투과 또는 소멸되는 특징을 갖는다. 따라서, 사석구조물의 성능을 정확히 파악하기 위해서는 제체의 외부뿐만 아니라 공극매체 내부에서의 파동장을 동시에 해석하여야 한다.

공극매체내에서의 파동장 해석을 위해서는 물과 토립자의 혼합영역에서의 지배방정식 및 경계조건들을 먼저 검토하여야 한다. 사실 상 공극매체를 하나의 연속체로 고려하여 일종의 Navier- Stokes 방정식(이하 N-S 방정식으로 칭함)을 수립하고자 하는 시도는 역사가 꽤 오래 되었으며 Biot (1956), De Josselin de Jong (1956), Pride et al.(1992), Verruijt(1994)가 관련연구를 수행한 바 있다. 공극매체에서의 파동장 해석을 위해서는 Sollitt and Cross (1972)가 N-S 방정식에 근거하여 사석방파제의 반사와 투과를 이론적으로 해석하기 위한 연구를 수행하였다. 그 이후 van Gent (1995) 역시 공극매체내에서의 N-S 방정식을 유도하였으며 VOF 방법을 이용하여 사석방파제에 적용한 바 있다. Cruz et al. (1997)은 공극매체내에서의 Boussinesq 방정식을 유도하고 이를 수중 잠제의 성능을 검토하기 위하여 적용하였으며, Lynett et al. (2000)은 역시 형태는 다르나 일종의 Boussinesq 방정식을 유도하여 고립파 입사에 대한 사석방파제의 성능을 검토한 바 있다. 한편, Sakakiyama and Kajima (1992)는 공극매체내의 N-S 방정식을 직접 해석하는 일종의 RANS 모델을 제시하였으며 현재의 CADMAS-SURF 프로그램(CDIT, 2001)의 기초를 제공하였다. 최근엔 Hu et al. (2012)이 dam-break 파와 공극체간의 상호작용에 대한 수치해석을 수행하였으며 국내에서는 이창훈·박설화 (2011)가 완경사방정식을 이용하여 투수층에서의 불규칙파 전파를 모의한 바 있다.

상기와 같이 공극매체내에서의 파동장에 대한 연구가 지속적으로 수행되어 왔음에도 불구하고 관련 문헌들을 검토하여 보면 N-S 방정식에 근거한 지배방정식을 구성하는 각 항들에서 공극률의 적용 형태가 다소 다름을 발견할 수 있다. 일반적인 물 영역에서의 파동장과는 달리 공극매체내에서는 토립자에서의 마찰력과 관성력에 따른 모멘텀 교환, 즉, 공극 항력이 전체의 파동장 결정에 있어서 매우 중요한 역할을 한다. 따라서, 지배방정식에서 각 항들의 비교에서 공극 항력의 상대적 크기가 잘못 설정될 경우 계산결과가 현격히 다르게 나타날 수 있다.

본 연구에서는 연속 공극매체 내에 미소부피의 검사체적을 설정하고 Reynolds 이송정리를 적용하여 질량과 모멘텀의 교환을 고려하는 전통적인 방법에 따라 N-S 방정식을 유도하고 이를 기존 문헌들에서의 식들과 비교하였다. 그리고 유도된 N-S 방정식 및 경계조건들을 이용하여 공극매체내의 Boussinesq 방정식을 유도하였으며 이차원 직립사석 방파제에 대하여 수치해석을 수행하고 이 결과를 기존의 수리실험결과들과 비교하였다.

2. 공극매체 Navier-Stokes 방정식의 유도 및 검토

Fig. 1에서 보이는 바와 같이 공극률(porosity)이 인 등방성 공극매체로 충진되어 있으며 체적과 질량이 각각 와 인 고정된 미소검사체적에서의 임의 스칼라 또는 벡터에 대하여 Reynolds 이송정리를 적용하면 다음과 같다.

(1)

여기서, 는 검사표면, 는 임의 시간에 에 존재하는 임의 시스템(질량, 모멘트, 에너지 등), 는 단위질량당 시스템(), 는 단위 외향법선벡터, 는 침투유속(Seepage velocity) 벡터를 각각 의미한다.

Eq. (1)에서 유체흐름의 연속방정식과 운동방정식은 각각 질량()과 모멘텀()의 이송을 나타낸다. 직각좌표계(연직상향 +)에서 비압축성 흐름의 경우에 이들 방정식을 지표(Index)형태로 표현하면

;        (2)

 ;

       (3)

와 같다. 여기서, 는 차례로 축을 의미한다. Eq. (3)의 우측 항들은 차례대로 에 가해지는 표면력(압력+점성마찰력), 자체 중력, 공극체 입자에 의한 내부마찰력(Forchheimer 식), 그리고 공극체 입자에 작용하는 관성력을 의미하며 각각 다음과 같이 나타낼 수 있다.

 (4)

 (5)

 (6)

 (7)

Eq. (4)에서 는 점성계수, (5)에서 는 중력가속도 성분 을 의미한다. Eq. (6)에서 , 는 Forchheimer 식에서의 공극율과 토립자 입경에 의존하는 유차원 파라미터들이다(단위는 각각 [], [] 임). Eq. (7)에서 은 가상질량계수를 의미하며 부가질량계수 와 의 관계가 있다. Eq. (6)에서 는 Darcy 유속으로서 침투유속 와는 다음과 같은 관계가 있다.

 (8)

Eqs. (4)~(7)을 Eq. (3)에 대입하면, 최종적으로 공극매체에 대한 운동방정식을

 (9)

와 같이 나타낼 수 있다. 여기서, , , 이며

 (10)

이다. Eq. (9)를 Darcy 유속 를 사용하여 다시 쓰면 다음 식과 같다.

 (11)

Eq. (11)에서 정상 층류 가정 하에 가속도항과 내부마찰력 항을 무시하고 공극마찰력에서 층류항만을 고려하면

 또는   (12)

이 되어 결국 Darcy의 법칙과 일치하게 된다.

공극마찰에 대한 Eq. (11)에서의 계수는

,     (13)

와 같다(단위는 각각 [], [] 임). 여기서, 는 고유투수성 (Intrinsic permeability), 는 공극 내 흐름 통로의 형상에 관계되는 무차원수를 의미하며 통상 0.6의 값을 갖는다. Fair and Hatch (1933)는 고유투수성을 다음과 같이 제시하였다.

 (14)

여기서, 5.0, 7.7이며 은 공극매체의 평균직경이다. 한편, Engelund (1953)는 Eq. (13)의 , 를 다음과 같이 제시하였다.

 (15)

 (16)

Eqs. (15) and (16)에서 = 780~1500 또는 그 이상, = 1.8~3.6 또는 그 이상이다.

운동방정식 (9)와 함께 기존의 타 연구들에서 사용된 식 들을 Table 1에 제시하였다. 각 식에서 유속은 침투유속 로 통일하였다. Sollitt and Cross (1972)의 식은 이류가속도 항이 빠져 있으나 Eq. (9)와 동일하며, Sakakiyama and Kajima (1992)의 식도 내부마찰항이 항력계수 로 표현된 것을 제외하면 등방성 매체에서 이므로 Eq. (9)와 일치한다. 그러나, Cruz et al. (1997)의 식은 거의 Eq. (9)와 유사한 형태를 보임에도 불구하고 좌변에서 계수 에 공극률이 곱해져 있다. 한편, van Gent (1995)의 식에서는 Eq. (9)와 비교하여 이류가속도항에서 차이가 발견되며 내부마찰항도 달리 표현되어 있음을 볼 수 있다.

Burcharth (1995)는 관성력항과 압력항의 관계에 대한 고찰에서 Eq. (9)의 형태보다는 Cruz et al. 식의 형태가 더 합당한 것으로 판단한 바 있다. 그러나, Cruz et al.의 식이 의미를 갖기 위해서는 Eqs. (4), (5) and (6)에서 가 아닌 가 되어야 한다. 이는 Reynolds 이송정리를 이용한 본 유도과정에서 알 수 있듯이 전혀 타당하지 않다. 다만, Cruz et al.에서의 계수 를 또 다른 계수 으로 치환하면 형태상으로는 Eq. (9)와 동일하게 된다. 이는 Cruz et al.의 식이 타당성을 갖기 위해서는 계수 이 Eq. (10)에 의하여 결정되는 이론적 계수가 아니고 일종의 실험적으로 결정되는 계수라야 한다. van Gent의 식에서도 계수 는 실험적 상수 에 의존하는 일종의 실험적 계수로 표현되어 있으며 가상질량계수 이 부가질량계수 와 의 관계로 인하여 1보다 클 것이라는 예상과는 달리 으로 오히려 1보다 작게 나타나 있다. Sakakiyama and Kajima (1992)도 역시 를 실험적으로 결정하여 사용하였다. 이와 같은 사실을 종합해볼 때, 공극매체 내부에서 토립자에 작용하는 관성력은 고립된 단일 입방체에 작용하는 관성력과는 매우 다른 작용 메커니즘을 갖는 것으로 보인다. 이에 대해서는 차후 수치해석에서 다시 논의하기로 한다.

3. 공극매체에서의 Boussinesq 방정식의 유도 및 적용

3.1 Boussinesq 방정식의 유도

Fig. 2의 이차원( 평면) 공극매체 파동장에 비압축성, 비점성의 가정을 도입하면 연속방정식과 운동방정식은 Eqs. (2) and (9)에서의 을 로 치환하여 Eqs. (17)~(19)와 같이 나타낼 수 있다.

 (17)

 (18)

 (19)

여기에 비회전 유체흐름의 가정을 도입하면,

 (20)

이다. Eqs. (17)~(20)을 풀기위한 경계조건은 다음과 같다.

- 운동학적 자유수면 경계조건

 on   (21)

- 동역학적 자유수면 경계조건

 on    (22)

- 바닥경계조건

 on       (23)

Boussinesq 방정식은 파동의 진폭 , 대표수심 그리고 파장 사이의 두개의 무차원 파라미터, 즉 비선형파라미터 와 천수파라미터 에 대하여 Eqs. (17)~(23)을 섭동시켜 전개할 수 있다. 비 공극매체 흐름에 대해서는 Nwogu (1993)와 Wei et al. (1995)이 속도포텐셜을 미지수로 하는 지배방정식을 이용하여 Boussinesq 방정식을 유도한 바 있다(Dingemans, 1997). 한편, Walkley (1999)는 직접 Eqs. (17)~(19)에 상응하는 비공극매체에서의 지배방정식을 이용하여 Peregrine (1972)의 식과 Nwogu (1993)의 확장형 Boussinesq 식을 각각 유도하였다. 전자는 적용범위가 천해파에 국한되는 반면, 후자는 수심변화에 따른 분산성의 재현정도를 심해파로 확장한 것이다.

본 연구에서는 Walkley (1999)의 방법에 따라 공극매체에 대한 Nwogu 타입의 Boussinesq 방정식을 유도하기로 한다. 먼저, 파라미터 , 와 를 이용하여 Eqs. (17)~(23)의 모든 시스템 변수들을 무차원화하면 다음과 같다.

, , , ,

, ,

Boussinesq 방정식의 스케일해석을 위하여 과 를 이용하여 , , , , 의 자릿수를 으로 조정하면 다음과 같다.

, , ,

,

상기 무차원수들을 Eqs. (17)~(23)에 각각 대입하면 다음의 Eqs. (24)~(30)을 얻을 수 있다.

 (24)

 (25)

 

(26)

 (27)

 (28)

 (29)

 (30)

연속방정식 (24)로 부터

 (31)

이며 Eq. (27)에 대입하면

 (32)

와 같다. 임의점 에서의 수평유속을 라고 정의하고 를 에 대하여 Taylor 전개하면

(33)

와 같으며, 에서 까지 적분하면 다음과 같다.

(34)

Eq. (34)를 Eq. (32)에 대입하면,

      

(35)

이며 이를 다시 에 대해서 미분하면, 

 (36)

와 같다. Eqs. (35)와 (36)을 Eq. (33)에 대입하면,

(37)

이며, Eq. (37)을 Eq. (31)에 대입하면

 (38)

을 얻을 수 있다.

Eqs. (37) and (38)을 연직모멘텀식 (26)에 대입하여 를 구한 다음, 동역학적 자유수면 경계조건인 Eq. (29)를 이용하여 에서 자유수면 까지 적분하면 압력의 연직분포를 다음과 같이 구할 수 있다.

 

(39)

Eqs. (37)~(39)를 수평모멘텀식 (25)에 대입하면 자릿수가 의 크기인 다음과 같은 첫 번째 Boussinesq 방정식인 운동방정식을 얻을 수 있다.

 (40)

연속방정식 (24)를 연직적분하고 운동학적 수면경계조건 (28)을 적용하면

 (41)

를 얻을 수 있으며 여기에 Eq. (37)을 대입하면 다음과 같은 두 번째 Boussinesq 방정식인 연속방정식을 얻을 수 있다.

 (42)

       

Eqs. (42) and (40)을 유차원으로 환원하고 로 정의하면 (), 공극매체에서의 2D Boussinesq 방정식은 다음과 같다.

(43)

  

 (44)

여기서, , , , 이다. 상기 식들을 Darcy 유속 으로 표현하고 Eq. (44)에 Eq. (11)의 내부마찰항을 부가하면 다음과 같다.

 (45)

   (46)

Eqs. (45) and (46)은 공극률이 이면 순수 유체흐름에 대한 Walkley (1999)와 결과와 정확히 일치한다(기존의 Dingemans (1997, p.537)의 결과와도 일치하나 Nwogu (1993)의 결과와는 Eq. (46)의 4, 5번째 항에서 약간의 차이가 발생함). 한편, Nwogu는 적정 값으로서 을 제안하였다.

Eqs. (45) and (46)을 삼차원으로 확장하면,

 (47)

 

 (48)

와 같다. 여기서, 는 에서의 수평유속 성분을 나타낸다. Wei et al. (1995)는 잔류항이 인 Eqs. (47) and (48)의 비선형성을 보다 확장하여 잔류항을 천수파라미터 가 주도하는 의 연속방정식과 운동방정식을 유도하였다.

Boussinesq 방정식에 대한 수치해석은 Wei et al. (1995, 1999), Nwogu and Demirbilek (2001)에 제시되어 있다. 수치해석의 안정성을 위하여 시간 및 공간차분의 증분을 일반적인 Courant 조건에 따라 결정하면 되나 본 확장형 방정식에서 모든 격자에서의 수심조건을 분산성의 심해측 유효한계인 수심/파장=0.5 이내로 유지하는 것 그리고 한 파장내에 최소 8개 이상의 격자를 두는 것 등에 유의하여야 한다. 특히, 불규칙파 입사조건을 채택하는 경우는 고주파수 성분에서의 해상도 및 수치해석의 안정성을 위하여 시간 및 공간증분, 전체 계산시간 등을 세심히 고려할 필요가 있다.

3.2 직립식 공극 방파제에 대한 적용

유도된 공극매체 Boussinesq 방정식의 적용성을 검토하기 위하여 직립식 공극 방파제의 반사율과 투과율을 각각 계산하여 기존의 수리실험결과와 비교하였다. 수치해석은 Eqs. (45) and (46)에 Wei et al. (1995)의 Predictor-Corrector 방법을 적용하였으며 기존의 Boussinesq 방정식 모델(Chun et al., 2006; Chun, 2007)을 수정하여 사용하였다. 계산에 사용된 수치수조는 Fig. 3에서 보이는 바와 같이 수조의 좌·우단에 파를 흡수하기 위한 스펀지층을 설치하였으며 반사율과 투과율을 기록하기 위한 수치파고계들을 공극방파제 전·후에 설치하였다. 아울러, 해석결과의 적절한 평가를 위하여 기존에 수립된 경계요소법 모델(Chun et al., 2007)도 함께 적용하였다. 본 경계요소법에서는 Eq. (11)에서 이류항을 무시하고 가상질량계수를 =0으로 하여 적용한 것이다.

수치해석의 입력조건은 Table 2에서 보이는 바와 같이 Sollitt and Cross (1976)와 Keulegan (1973)의 수리실험들과 동일하게 취하였으며 계산결과를 이들의 수리실험결과들과 비교하였다. 본 수리실험들은 고정된 주기에 파고를 증가시켜가며 파고계를 이용하여 반사율과 투과율을 측정하였다. 실제 본 수리실험들은 주기와 파고를 달리하는 다양한 조건에서 수행되었으나 본 연구에서는 현 수치해석기법의 성능을 일차적으로 검토할 목적으로 Table 2와 같이 각각 한 주기들만을 선별하여 수치해석을 수행하였다.

수치해석에서 반사율은 Park et al. (1992)의 삼점법을 적용하여 계산하였으며 투과율은 zero-up-crossing 기법을 적용하여 계산하였다. 파고계 계측시간은 조파기 - 수조우단의 거리를 군속도로 3 왕복한 이후부터 6 왕복할 때까지로 하였다. 계산조건으로서는 = 으로 하였으며 시간증분은 Courant 조건을 적용하여 와 같이 취하였다. 여기서, 는 선형파 파속을 의미한다. Eq. (15)와 Eq. (16)에서의 내부마찰계수는 각각 = 1150 , = 2.7 로 취하였다.

본 수치해석에서는 Eq. (10)에서의 의 값에 대한 수치실험을 수행하였다. 즉, 쇄석에 대한 부가질량계수 를 동일 직경을 갖는 구체의 이론적 값 에서 관성력을 무시했을 때의 , 그리고 에 이르기 까지 다양한 값을 시도하여 계산결과를 수리실험결과와 비교하였다. Figs. 4 and 5는 Sollitt and Cross (1976)에 대한 계산결과를 수리실험결과와 비교한 것이다. 반사율과 투과율 공히 일 때 수리실험결과와 가장 잘 일치함을 볼 수 있다. 이는 에서 부가질량계수가 이 됨을 의미한다. 이와 같이 음의 값을 갖는 부가질량계수가 물리적으로 전혀 불가능한 것은 아니나 본 문제에서와 같은 극단적인 값 이 물리적인 타당성을 갖는다고 보기는 매우 어려울 것으로 보인다. 대신 공극매체 흐름에서 쇄석부분에 의한 관성저항은 형성되지 않으며 형성된다 하더라도 거의 무시될 수 있는 것으로 간주하는 것이 타당할 것으로 판단된다. 한편, 경계요소법의 결과는 수리실험결과와 비교적 잘 일치하는 것으로 나타나 있다. 본 경계요소법의 결과도 의 조건하에서 계산된 것이다(Sollitt and Cross (1972)도 을 특정하기가 어려운 점을 고려하여 로 하여 수치해석을 수행한 바 있다). 또한, Fig. 4의 반사율에서 에 대한 Boussinesq 모델의 결과와 경계요소법 결과가 전 의 범위에서 거의 일정한 차이를 보이고 있으며 변화 성향이 거의 동일함을 볼 수 있다. 이는 경계요소법에서 이류항이 제외되었기 때문인 것으로 판단된다.

Fig. 6은 의 조건에서 파랑의 1 주기 내에서 네 개의 동일 시간간격마다 계측한 전 수조연장(150 m)에서의 스냅 파형이다. 공극방파제 전면에는 부분 중복파의 형성에 따른 파형의 시간적 변화를 볼 수 있으나 후면에는 일정 진폭의 투과파가 형성되어 있음을 볼 수 있다.

Keulegan의 수리실험 조건에 대한 계산결과를 수리실험결과와 함께 Figs. 7~10에 도시하였다. 상기 Sollitt의 실험조건에 대한 결과에서와 마찬가지로 의 조건에서 일치정도가 가장 양호하게 나타났다. 그러나, 전반적으로 투과율에서는 계산치와 실험치가 비교적 잘 일치하는 반면 반사율에서는 상당 정도의 편차가 존재함을 볼 수 있다. 반사율에서 대체적으로 Boussinesq 모델결과와 경계요소법의 결과가 유사한 것으로 미루어보아 본 계산치와 실험치와의 차이는 반사율의 측정방법의 차이에 의하여 발생한 것이거나 아니면 공극매체에서의 토립자에 의한 내부마찰력의 차이에 기인한 것으로 판단된다.

Keulegan의 수리실험에서는 1개의 파고계를 이동시켜가며 얻은 시그널로 반사율을 구하였지만 본 Boussinesq 수치해석에서는 3개의 독립된 파고계의 시그널를 이용하는 삼점법을 사용하였다. 내부마찰력에 대해서는 Eq. (11)에서 알 수 있듯이 수립자의 속도에 의존하기 때문에 수치해석의 정확도는 수립자 궤도유속(orbital velocity)의 연직분포를 얼마나 잘 재현하는가에 달려있다. 그러나, 본 수치해석에서 사용한 확장형 Boussinesq 모델은 주로 심해역에서의 분산성의 재현 정도를 개선하기 위한 것으로서 수립자의 궤도유속을 잘 재현한다고 보기는 어렵다. 일단, Eqs. (37) and (38)에서 수립자 유속의 오차크기를 로 작게 제시하고는 있으나 분산성과 비선형성을 달리하는 여러 파랑조건들에 대한 별도의 적극적인 검증이 필요할 것으로 판단된다. 한편, 수심을 여러 층으로 나누어 수립자 궤도유속 계산의 해상도를 높이기 위한 기법들이 최근 연구된 바 있다(Lynett and Liu, 2004; Hsiao et al., 2005). 이들 기법들을 본 수치해석에 사용함으로써 계산의 정확도를 제고할 수 있을 것이다. 아울러, 수리실험 역시 매우 오래전에 수행되었기 때문에 최소한 몇 개 조건에 대해서는 반사율 및 투과율의 확인실험이 수행될 필요가 있을 것으로 판단된다.

4. 결 론

공극매체 흐름에 대한 Reynolds 이송정리를 적용하여 공극매체에서의 Navier-Stokes 방정식을 유도하였으며 이를 기존의 유사연구 결과들과 비교하였다. 또한, 유도된 N-S 방정식을 이용하여 공극매체 내외에서 파동장의 비선형성과 분산성을 적절히 재현하기 위한 확장형 Boussinesq 방정식을 유도하였다. 그리고 이들 방정식의 유효성을 검증하기 위하여 공극방파제의 반사율 및 투과율에 대한 수치해석을 수행하여 그 결과를 기존 수리 실험들과의 비교하였다. 세부적인 연구결과는 다음과 같다.

(1)수치실험을 수행한 결과, 공극방파제의 반사율과 투과율은 토립자의 가상질량계수에 민감하게 반응하는 것으로 나타났으며, 가상질량계수를 영으로 처리했을 경우에 수리실험결과와의 일치정도가 가장 양호한 것으로 나타났다. 이는 공극매체 파동장에서 토립자에 의한 저항은 주로 내부마찰로 나타나며 관성저항은 매우 미미하게 형성됨을 의미한다.

(2)수치해석결과에서 투과율은 수리실험결과와의 양호한 일치를 보이는 반면 반사율에서는 다소 차이가 발생하였다. Boussnesq 방정식 모델의 결과가 경계요소법의 적용결과와 전반적으로 유사하게 나타난 것으로 보아 수리실험결과와의 차이는 반사율 계산방식의 차이에 근거하거나 공극매체내부에서 수립자의 수평 및 연직 궤도유속의 재현정도가 다소 미흡하기 때문인 것으로 판단된다. 이를 개선하기 위한 추가적인 연구 및 확인 수리실험이 요구된다.