1. 서 론
지구온난화에 따른 해상 기상조건 악화와 더불어 전 지구적 무역량의 증가로 인한 대형 선박의 입출항 증가에 따라 해상 오염물질 유출사고가 빈번하게 발생하고 있다(Xie et al., 2017). 연안에서의 스칼라 이송에 대한 이해는 오염물질 유출(Rivord et al., 2014), 유류유출(Jung and Son, 2018), 녹조(Rowe et al., 2016) 등 연안에서 발생 가능한 해양오염 재해를 해석하고 방지하기 위해 필수적이다. 염료 등 추적자를 이용한 연안에서의 직접적인 관측은 실험 규모의 제약 등 다양한 문제를 수반하여 제한적으로 수행되고 있기 때문에 일반적으로 수치모형이 폭넓게 활용되고 있다.
비선형 천수방정식 수치모형은 높은 계산 효율성과 안정성으로 조석(Rahman et al., 2020), 폭풍해일(Lee et al., 2019), 지진해일(Son et al., 2011) 등 다양한 장주기파의 해석에 주로 사용되고 있다. 또한, 비분산성 및 정수압 조건에서의 넓은 영역의 흐름과 스칼라 이송 예측을 위해 수심평균된 이송확산방정식과 결합하여 일반적으로 사용된다(Kim and Lynett, 2013). 그러나 비선형 천수방정식 수치모형은 ‘완벽한’ 장파의 모의만 가능하며, 파랑의 분산성을 고려할 수 없기 때문에 연안에서 파랑의 불규칙한 특성을 반영하는 데 한계가 있다.
Boussinesq 파랑모형은 파랑의 분산성을 고려할 수 있으며 장파뿐 아니라 단파 및 불규칙파의 모의가 가능하기 때문에 다양한 흐름의 해석에 적용할 수 있어 많은 연구가 진행되고 있다. 더불어 파랑-흐름간 상호작용을 고려할 수 있으며, 쇄파로 인해 발생하는 연안류의 해석이 가능하다는 장점이 있다(Son and Lynett, 2014; Choi and Seo, 2018). 그러나 파랑의 분산성을 고려하기 위한 고차항의 계산은 연산량을 증가시킴과 동시에 수치모형의 안정성을 저하한다. 수치모형의 안정성을 증가시키기 위한 음해법의 사용이 추천되나, 연산량을 더욱 가중시키는 단점을 야기한다. 따라서, Boussinesq 파랑모형의 사용은 빠른 계산을 위한 하드웨어 구축에 많은 비용을 수반하며, 수백 개의 CPU 코어를 활용한 병렬연산, GPU를 이용한 연산 등 고성능 컴퓨팅시스템을 요구하게 된다(Erduran et al., 2005; Tavakkol and Lynett, 2017).
GPU는 수치모형의 연산을 가속화할 수 있는 많은 계산자원 중 경제적으로 합리적인 대안이기 때문에 이를 활용한 수치모형의 개발에 관한 많은 연구가 수행 중에 있다(Gandham et al., 2015; Tavakkol and Lynett, 2017; Yuan et al., 2020). Gandham et al. (2015)은 GPU 가속화된 불연속 갤러킨 기법을 이용한 천수방정식 수치모형을 개발하였으며, Yuan et al. (2020)은 기개발된 CPU 기반의 Boussinesq 파랑모형인 FUNWAVE-Total Variation Diminishing (TVD)을 다중 GPU 코어를 이용하여 가속화한 FUNWAVE-GPU를 개발하였다. 마찬가지로, Tavakkol and Lynett (2017)은 GPU 가속화를 활용하여 Boussinesq 방정식을 풀이하는 파랑모형인 Celeris Advent를 개발하였다. Celeris Advent는 수치적인 안정성을 증가시키기 위해 Minmod 제한자를 사용하였으며, Graphical User Interface (GUI)를 통해 수치모의 진행 중에 실시간으로 결과를 가시화할 수 있을 뿐만 아니라 마우스 입력을 이용하여 수위, 수심지형을 바꿀 수 있는 상호작용 시스템을 갖춘 최초의 수치모형이다.
본 연구에서는 GPU 가속화를 활용한 Boussinesq 수치모형인 Celeris Advent에 이송확산방정식을 결합하여 흐름 및 스칼라 이송-확산 예측모형을 개발하였다. 더불어 1차원 및 2차원 공간에서의 스칼라 이송의 해석해와 실험결과를 수치모의 결과와 비교함으로써 개발모형의 정확도를 검증하였다.
2. 지배방정식
Celeris Advent의 지배방정식은 Madsen and Sørensen (1992)의 확장형 Boussinesq 방정식을 사용하고 있으며, Tavakkol and Lynett (2017)은 이를 보존형 방정식으로 변환하여 수치적으로 이산화하였다. 본 연구에서는 3차원의 이송확산방정식을 수심에 대해 적분한 수심평균된 이송확산방정식을 지배방정식 내에 추가하였으며 이에 따른 지배방정식은 Eq. (1)과 같다.
여기서, h는 전수심, t는 시간, (P, Q)는 각각 x, y 방향의 수심적분된 질량 플럭스, c는 스칼라의 농도, g는 중력가속도, (u ¯ v ¯
3. 수치해석 기법
3.1 Hybrid Finite Volume-Finite Difference Method
유한체적법은 플럭스 재구성와 제한자의 사용에 따라 수치모형의 안정성을 증가시키는 장점을 가지고 있다. 천수방정식과 달리 Boussinesq 방정식은 분산항의 고려에 따라 유한체적법을 통해 쉽게 수치적으로 이산화되지 않는다. Erduran et al. (2005)은 Boussinesq 방정식의 수치이산화를 위해 유한체적법과 유한차분법을 동시에 사용하는 Hybrid finite volume-finite difference method를 제안하였다. 제안된 Hybrid scheme은 Boussinesq 방정식을 풀이함에 따른 높은 해석 정확성과 함께 유한체적법을 사용함에 따라 수치모형의 안정성을 증가시키는 장점이 있어, 많은 선행연구에서 이를 적용하여 수치모형을 개발하였다(Shi et al., 2012; Kim and Lynett, 2013; Tavakkok and Lynett, 2017). Erduran et al. (2005)에서 제안한 방법과 같이 지배방정식의 이송항에 대해서는 유한체적법을, 분산항, 확산항 등 나머지 항에 대해서는 유한차분법을 바탕으로 각각 이산화하였다. Kurganov and Petrova (2007)가 제안한 2차 Well-balanced positivity preserving central-upwind 수치기법을 이용하여 격자점 경계에서의 플럭스Hi는 Eq. (2)와 같이 계산한다.
여기서, 아래첨자i는 격자점의 경계지점을 의미하며, (nx, ny)는 격자점 경계에서의 x, y 방향의 법선벡터를 의미한다. a i ± ∂ F ∂ U
3.2 시간적분법
본 연구에서는 Tavakkol and Lynett (2017)에서 적용한 것과 동일하게 시간미분항의 계산을 위해 예측자-수정자 방법을 사용하였으며, 예측자 기법은 3차정확도의 Adams-Bashforth 기법을 사용하였으며 Eqs. (3)~(6)과 같다.
여기서, d는 수심을 나타내며 아래첨자 x, y방향으로의 편미분을 의미한다. 다음 계산시간단계에서의 전수심 hn+1는 Eqs. (3), (7)을 통해 직접 계산할 수 있지만 운동량 Pn+1, Qn+1은 음함수 형태의 Eqs. (19), (20)을 풀이하여 계산할 수 있다.
여기서,
4. 수치모의 결과 및 검증
본 연구에서 개발한 이송-확산 예측모형의 정확도를 검토하기 위해 1, 2차원 공간에서의 흐름 및 파랑조건에 대해 수치모의을 수행하였다. 첫 번째 실험에서는 1차원 일정수심의 수로 내 등류 흐름에서스칼라의 이송을 모의하였으며, 두 번째 실험에서는 원형 천퇴 지형의 후류에서 발생하는 와류에 의한 스칼라 이송을 모의하였다. 마지막 실험에서는 3차원의 복잡한 지형을 통과하는 쇄파 고립파의 전파에 의한 스칼라 이송을 모의하였다.
4.1 1차원 등류에 의한 스칼라 이송 모의
길이 16 m의 마찰이 없는 수로 내 유속 1 m/s를 가정하였으며, 스칼라의 초기조건은 다음과 같다.
격자크기∆x=0.01 m를이용하여 1,601개의 격자점을 구성하였으며, ∆t=0.001 sec를 이용하여 총 모의시간 10초 동안의 흐름에 의한 스칼라 이송을 계산하였다. 이송항의 수치적 계산시 발생하는 수치확산을 확인하기 위해 확산항은 고려하지 않았다.
Fig. 1은 이송-확산 예측모형을 이용한 계산결과와 해석해와의 비교로, 계산된 결과는 해석해와 매우 정확하게 일치하였다. 또한, 스칼라가 이송되는 동안 수치확산은 거의 발생하지 않았음을 알 수 있다.
4.2 2차원 원형 천퇴 지형의 후류에 의한 스칼라 이송 모의
Lloyd and Stansby (1997)는 개수로에서 등류 흐름이 원형 천퇴 지형을 통과하며 발생하는 후류에서의 와류방출 변화를 모의하기 위한 수리실험을 수행하였다. 상류에서 다양한 유속조건, 개수로 내 수심조건 및 천퇴 지형의 형상을 사용하여 실험을 수행하였으며, 메틸렌블루 염료를 이용하여 후류에서의 와류방출에 따른 스칼라 이송 또한 관측하였다.
본 연구에서는 Lloyd and Stansby (1997)의 수리실험 조건 중 SB4_02를 수치적으로 재현하였으며, 수치모의에 사용된 원형 천퇴가 있는 수심 지형은 Fig. 2에 제시되어 있다. 길이 9.75 m, 폭 1.52 m의 수로의 유입구로부터 5.0 m 하류 부근에 반경 0.75 m의 원형 천퇴가 설치되었다. 천퇴를 제외한 영역의 수심은 0.054 m이다. 천퇴의 높이는 0.049 m이며, 천퇴의 경사각은 8.0°이고, 천퇴의 중심부에서의 수심은 0.005 m로 천퇴의 높이와 수로의 수심의 비율은 1.10이다. 수로 내 상류의 유속조건은 0.115 m/s로 일정하며, 하상에서의 마찰을 고려하기 위해 매닝 조도계수를 0.02로 적용하였다.
격자크기는 ∆x=∆y=0.015 m를 사용하였으며, 스칼라는 천퇴의 전면 사면에서 주입되는 것으로 가정하였다. 수치모형을 통해 수리실험이 잘 재현되었는지 확인하기 위해 천퇴의 중심부에서 하류 방향으로 1.02 m 떨어진 지점에서 관측된 유속을 비교하였다. 관측점 A는 y축 방향으로 수로의 중심선에 위치하며, 관측점 B는 중심선으로부터 0.27 m 떨어진 지점에 위치한다.
Fig. 3은 관측점 A, B에서의 개발모형을 이용한 계산결과와 관측치를 비교한 결과를 나타낸다. 유속 v의 경우 두 관측점에서 관측치와 최대값과 상이 잘 일치함을 확인되었다. 그러나, 유속 u는 불규칙한 패턴을 보인 실험결과와 달리 계산결과는 규칙적인 패턴을 나타냈다. 또한, 계산결과는 실험결과와 비교했을 때 전체적으로 과소추정되었으며 관측점 A에서 차이가 더욱 크게 발생하였다. 이러한 결과는 수심적분형 수치모형을 사용하여 본 수리실험을 재현한 선행연구(Lloyd and Stansby, 1997; Zhang et al., 2016)에서도 공통적으로 발견되었으며, 이는 수심적분형 수치모형이 복잡한 3차원 지형조건에 의해 발생하는 난류를 충분히 반영하지 못하여 발생하는 것으로 확인되었다. 그러나 선행연구에서 수행한 수치모의 결과와 비교했을 때 유속의 양상은 잘 일치하는 것을 확인하였다.
Fig. 4(b)는 개발모형을 통해 계산한 순간적인 스칼라 농도 분포로 Lloyd and Stansby (1997)의 수리실험에서 관측한 스칼라 농도 분포인 Fig. 4(a)와 유사한 양상을 보여주고 있다.
Fig. 4
Instantaneous Scalar Distribution. (a) by loyd and Stansby (1997) Experiment and (b) Using Developed Model
4.3 불규칙한 수심 지형을 따라 전파하는 고립파에 의한 스칼라 이송 모의
Lynett et al. (2019)은 고립파가 불규칙한 수심 지형을 따라 전파하는 동안 발생하는 난류, 쇄파, 처오름 등 복잡한 흐름특성을 분석하고자 수리실험을 수행하였다. 고립파가 전파하는 동안의 수위를 측정함과 동시에 원뿔형 섬 후면의 고정된 지점에서 염료를 방출하여 염료의 이송 및 확산과정을 관측하였다.
본 연구에서는 Lynett et al. (2019)의 개수로에서의 고립파 전파 수리실험을 개발모형을 통해 수치적으로 재현하였다. Fig. 5는 수치모의에 사용된 Lynett et al. (2019)의 수리실험에서 사용된 수심 지형과 수위 관측지점의 위치를 나타낸다. 9개의 관측지점 중 y축의 중심에 위치한 4개의 관측지점만 사용하였으며, 4개의 관측지점은 x축에 대해 각각 7.5, 13, 21, 25 m에 위치하고 있다.
x = 1 m에서 수심의 절반에 해당하는 0.39 m 파고의 고립파를 발생시켰으며, x 방향의 모의영역의 양측 경계에 스펀지층을 설치하여 경계면에서의 반사파의 생성을 방지하였다. 격자크기는∆x=∆y=0.1 m를 사용하였으며 바닥 마찰을 고려하기 위해 바닥 마찰계수는 0.0025를 사용하였다. Celeris Advent는 쇄파항을 고려하지 않기 때문에 파랑 쇄파에 의한 운동량 소실과 유사한 수치적 소실을 만들고자 Minmod 제한자를 사용하였으며(Tavakkol and Lynett, 2017), 이때 계수는 1.7을 사용하였다. 마지막으로, 스칼라는 원뿔형 섬의 후면에서 고립파가 원뿔형 섬의 정점을 통과한 직후부터 주입되는 것으로 가정하였다.
Fig. 6은 고립파가 전파하는 동안 관측점 1, 2, 3, 4에서 개발모형에 의해 계산된 수위와 실험결과의 비교를 나타낸다. 쇄파가 발생하기 전에 위치한 관측점 1, 2에서의 수치모의 결과는 실험결과와 잘 일치하는 것을 확인할 수 있다. 쇄파 발생지점보다 뒤에 위치한 관측점 3, 4에서는 쇄파로 인한 운동량 소실과 이에 따른 수위 감소를 충분히 반영하지 못하여 실험결과보다 크게 계산되었다. 최대 수위 값에서는 약간의 오차가 발생하였지만 전반적으로 Lynett et al. (2019)의 수리실험은 수치적으로 적절히 재현되었음을 확인하였다.
Fig. 7은 고립파의 전파 동안의 관측된 시간별 염료 분포를 나타낸다. 6.2초일 때 고립파는 원뿔형 섬의 정점을 통과하며, 쇄파가 발생하기 시작한다. 고립파는 아직 염료가 주입되는 지점에 도달하지 못했기 때문에 염료는 이송되지 않고 정지되어 있다. 8초일 때 고립파는 염료가 주입되는 지점을 통과하며 염료의 이동이 시작된다. 이후 경사면을 따라 처오름이 발생한 뒤 다시 내려오며 초기 고립파의 진행방향과 반대인 흐름이 생성된다. 20초일 때 초기에 생성된 고립파로 인한 정방향의 흐름과 사면을 타고 내려온 역방향의 흐름이 염료가 주입되는 지점의 인근으로 모이게 되어 염료가 집중된다. 각 방향의 흐름에 따라 염료는 분리되어 이송되며 27.4초일 때의 염료의 분포는 Fig. 7(d)와 같다. Fig. 8은 수치모의를 통해 계산한 시간별 스칼라 분포를 나타낸다. 전반적으로 복잡한 수심지형에서의 고립파 전파는 수치적으로 잘 재현되었으며, 스칼라 이송의 수치모의 결과는 실험결과와 유사한 양상을 보여주었다.
5. 결 론
본 연구에서는 GPU 가속화를 활용한 Boussinesq 수치모형인 Celeris Advent를 이용하여 확장형 Boussinesq 방정식과 수심평균된 이송확산방정식을 풀이하는 스칼라 이송 예측모형을 제시하였다. 해석해 및 실험결과를 이용한 세 차례의 수치모의 결과 비교 및 검증은 제시된 모형이 다양한 흐름 특성과 복잡한 지형에 따른 스칼라의 이송을 적절히 모의할 수 있음을 확인하였다. 1차원 공간에서의 등류 흐름에 의한 스칼라 이송에 관한 수치모의 결과는 수치확산이 거의 발견되지 않았으며, 해석해와 매우 정확하게 일치하였다. 또한 개수로에서의 흐름 및 파랑에 의한 스칼라 이송에 관한 수치모의 결과는 실험결과와 잘 일치하는 것을 확인하였다.
하지만, 개발된 수치모형은 쇄파에 의한 에너지 소실 및 난류혼합을 직접적으로 고려하지 못하는 한계점도 지니고 있다. 따라서, 쇄파에 의한 영향을 능동적으로 고려하기 위해 Kennedy et al. (2000)과 같은 Eddy viscosity 형태의 쇄파 방법의 적용이 필요할 것으로 보인다. 실제 연안에서의 오염사고를 모의하기 위해서는 연안의 흐름조건을 재현하기 위한 내부조파와 스폰지 경계층 경계처리가 향후 요구되며, 증발 등 오염물질의 물리화학적 특성에 대한 추가적 고려도 보완이 필요한 부분이다.
PDF Links
PubReader
ePub Link
Full text via DOI
Download Citation
Print
