1. 서 론
움직이는 탱크 내부에 부분적으로 적재된 유체가 일으키는 유동현상인 슬로싱(Sloshing), 해수가 선박의 갑판을 때리는 현상인 그린워터(Green-water) 등과 이에 동반하는 유체 충격하중 예측에 대한 연구가 활발히 진행되어 왔다. 최근에는 LNGC(Liquefied natural gas carrier), LNG-FPSO(Liquefied natural gasfloating production storage offloading), FSRU(Floating storage and regasification unit)와 같은 신개념의 선박들을 포함한 모든 종류의 선박이 대형화 되고 있어 선박의 구조적 안전성 및 운동성능 확보를 위한 슬로싱 현상의 정확한 해석에 관한 필요성이 강조되고 있다. 이뿐만 아니라, 2011년 동일본 대지진시 발생하여 큰 피해를 야기했던 지진해일의 경우처럼 재난방재를 위해서도 유체충격하중의 정확한 예측 시뮬레이션 기술의 필요성이 증대되고 있다.
일반적으로, 자유표면을 포함하는 물리현상은 비선형성이 강하기 때문에, 수치적 해석이 매우 까다롭고 복잡하다. 오일러리안 프레임(Eulerian frame)을 기반으로 자유표면 근처에 라그란지안(Lagrangian) 입자나 세그먼트(Segment)를 배치하여 그 거동을 풀거나 이동격자를 이용하는 방법들은 전도(Turnover) 등의 현상이 일어나지 않는 비교적 완만한(Mild)한 단상류(Single phase flow) 해석에 주로 이용되어져 왔으나, 쇄파(Wave breaking)나 기체-액체의 비선형성이 강하게 상호작용하는 경우에는 적합하지 않다. 이외에 대표적인 격자 기반의 오일러 접근법으로는 VOF (Volume of fluid)법(Hirt and Nichols, 1981), Level-set법(Sussman et al., 1994), 밀도함수법(Miyata and Park, 1995) 및 이를 개량한 방법들이 다수 존재한다. 오일러리안 프레임 기법에서 나타는 지배방정식의 대류항 처리기법은 해의 정확성과 계산 안정성에 큰 영향을 끼친다. 해의 정도 향상과 계산의 안정화를 위해서 고차의 스킴과 상류화(Upwinding) 기법을 사용하는 것이 일반적이나, 상류화기법을 사용하면 수치확산의 영향으로 복잡한 자유표면의 비선형 운동을 장시간 계산할 경우, 자유수면의 형상을 정확하게 묘사하기 힘들며, 질량이 보존되지 않는 등의 문제가 있다. 반면, 유동전체를 라그란지안 프레임에서 처리하는 대표적인 방법으로는 SPH(Smoothed particle hydro-dynamics)법(Monaghan, 1988), MPS(Moving particle simulation)법(Koshizuka and Oka, 1996)등이 있다. 이런 방법들은 대류항이 존재하지 않기 때문에 수치적 확산이 나타나지 않는다는 점과, 복잡한 유동도 비교적 간단히 처리 가능한 장점이 있다. 최근, Lee et al.(2010)은 MPS를 개량하여 MPS에서 고질적으로 나타나던 압력진동 문제를 해결하고, 정도 높은 결과를 보였으나, 국소적인 계산정도 향상과 3차원 문제의 계산시간 단축이 과제로 남아 있다.
각 계산점에서의 물리량과 이의 공간구배(미분)도 이용하는 CIP(Cubic interpolated propagation or constrained interpolation profile)법(Takewaki and Yabe, 1987)은 수치확산이 적다는 점과 인접 2점만의 정보로 3차 정도의 정확도를 갖는다는 등의 장점이 있다. Yabe et al.(1991)은 지배방정식의 모든 대류항에 CIP법을 적용하고, 압축성 및 비압축성 유동을 동일한 지배방정식을 이용하여 해석하는, 다상류문제(Multi-phase flow problem)의 해법으로서 CCUP(CIP-combined unified procedure)법을 제안하였다. CIP법의 개량형으로는 3차 보간함수의 근사에 유리함수를 이용하는 RCIP(Rational function CIP)법(Xiao et al., 1996)과 보존형 지배방정식을 사용하는 CIP-CSL(Conserved aemi-lagrangian)법(Xiao and Yabe, 2001)과 이의 개량형인 CIP-CSL2, CIP-CSL3, CIP-CSL4법 등이 있다.
한편, 오일러리안 프레임의 다상류해석을 수행하는 경우, 계면 포착(Interface capturing)을 위해서 밀도함수 등으로 정의한 자유표면의 수송방정식을 푸는 방법이 일반적이다. 이 때, 자유표면의 수송방정식의 대류항에 사용되는 CIP기반의 방법들로는 (R)CIP법, Tangent함수 등의 변환함수(Digitizer)법이 있으며, Xiao et al.(2005)에 의해서 제안된 보존형의 VOF계 기법의 하나인 THINC(Tangent of hyperbola for interface capturing)법은 다차원에의 확장이 용이하고, 수치 확산을 크게 줄일 수 있다고 알려져 있다. 최근에 Kim et al.(2008)은 기존 CIP법에 비해서 안정적이고, 계산량이 적은 USCIP(Unsplit semi-lagrangian, CIP)법을 제안하였다. Park et al.(2011)은 비압축성을 가정한 Navier-Stokes 방정식을 이용한 RCIP법 기반의 CCUP알고리즘과 VOF-PLIC(Piece-wise linear interface calculation)과 비슷하며 코딩이 보다 단순한 THINC-WLIC(Weighted line interface calculation)기법(Yokoi, 2007)을 채용하여, 2차원 직사각 탱크내의 슬로싱 문제를 해석하여, PNU-MPS(Pusan national university-MPS)법(Lee et al, 2010) 및 실험 결과와의 비교를 통하여 타당성을 검증하였다.
본 연구에서는, Park et al.(2011)에 의하여 개발된 CCUP법 기반의 계산수법을 3차원으로 확장하여, 자유표면의 거동, 기체 및 액체의 유동과 이에 동반하여 발생하는 유체 충격하중에 관한 문제에 적용하였다. 비압축성을 가정한 Navier-Stokes 방정식을 지배방정식으로, 대류항에는 USCIP법을 사용하였다. 계면 포착을 위해서는 THINC-WLIC기법을 사용하였으며, 압력 연성에는 MAC( Marker-and-cell)형의 알고리즘인 다단계(Fractional step)법을 채택하였고, 압력의 포아송(Poisson) 방정식은 SOR(Successive over relaxation)법을 이용하여 구했다. 개발된 수치해석기법의 검증을 위해 직사각 탱크내의 슬로싱 문제와 붕괴파에 의해 사각구조물에 작용하는 충격하중을 시뮬레이션하여 실험 결과와의 비교를 수행하였다.
2. 지배방정식
원래 CCUP법은 압축성 및 비압축성 유체가 단상 혹은 다상으로 존재하는 경우의 유동을 하나의 동일한 지배방정식을 사용하여 푸는 방법이지만, 본 연구에서는 코딩의 간결성을 위하여 압축성에 의한 영향은 미미하다는 가정하에 비압축성 Navier-stokes방정식을 지배방정식으로 사용하였으며, 이 때 기체-액체 경계면의 포착에는 식 (3)과 같은 밀도함수의 수송방정식을 사용하였다(Hu, 2004).
단, ρ=
, 𝜇=
, m=
, 𝜇=, m=
여기서, u는 속도, t는 시간, ρ는 밀도, p는 압력, 𝜇는 점성계수, Sij는 점성응력, fi는 중력 등의 외력이다.
식 (3)의 밀도함수 ϕ는 기체와 액체의 경우 각각 0.0과 1.0으로 설정하였으며, 식 (2)에 나타나는 각 격자에서의 평균 밀도 및 점성계수는 각 상의 대표값과 밀도함수의 값을 이용하여 구한다. Navier-stokes 방정식과 밀도함수의 수송방정식의 대류항에는 각각 USCIP법과 THINC-WLIC법을 사용하였다. 전체 계산 알고리즘은 Fractional step법으로, 상세한 내용은 다음 장에 서술한다.
3. 수치해석기법
3.1 (US)CIP법
CIP법의 기본적인 아이디어는 임의의 물리량 f의 분포가 유속 u에 의하여 이동하는 대류문제를 해석할 때, 이의 공간구배 g도 이용하여 분포를 구성하는 방법이다. 식 (4)와 같은 1차원 대류문제를 등간격 격자계에서 푸는 경우, Fig. 1(a)와 같이 원형 그대로 이동하여야 하는 물리량의 분포는, 선형내삽으로 인하여 Fig. 1(b)의 대시선(Dashed line)과 같은 형태(Fig. 1(c)에서 점선)로 변형되어 인식되어 진다. 반면, 식 (5)와 같이 물리량의 공간구배(미분)도 동시에 이용하면, Fig. 1(c)의 대시선과 같은 원형과 유사한 분포를 얻을 수 있다.
일반적으로 CIP법에서는 두 개의 계산점(xi와 xi-1)사이의 프로파일을 3차 다항식으로 근사하며, xi의 위치에서 얻어진 n시간 스텝에서의 프로파일인 Fnn을 −uΔt만큼 이동시킨, 
=F(x−uΔt), 
=dF(x−uΔt)/dx가 다음 시간 스텝(n+1)에서의 값이 된다. 격자의 사이즈를 1이라고 할 때, 이전 위치 p에서 추정된 물리량 f의 값은 식 (6)과 같이 나타낼 수 있다.
=F(x−uΔt), =dF(x−uΔt)/dx가 다음 시간 스텝(n+1)에서의 값이 된다. 격자의 사이즈를 1이라고 할 때, 이전 위치 p에서 추정된 물리량 f의 값은 식 (6)과 같이 나타낼 수 있다.
여기서, 계수 c0와 c1은 아래와 같이 두 격자점에서의 물리량 fi와 구배 gi로 나타내어진다.
2차원 문제의 경우, Yabe et al.(1991)은 아래와 같은 수식을 사용하였다.
상기식의 10개의 계수는 격자 모서리 점에서의 4개의 물리량과 6개의 공간구배를 이용해 구한다.
상술한 CIP법 뿐만 아니라, 2차 이상의 다항식을 근사식으로 이용하는 경우, 수치적인 진동이 발생할 수 있다. 최근에 Kim et al.(2008)은 기존 CIP법에 비해서 안정적이고, 계산량이 적은 USCIP법을 제안하였다. USCIP법에서는 모든 구배정보를 이용하고, 여기에 2개의 항을 추가함으로서 기지의 변수와 다항식의 개수의 숫자를 일치시킨다(식 (9)).
2차원 및 3차원 USCIP법의 계수는 참고문헌(Kim et al, 2008)에서 찾을 수 있다.
3.2 Fractional Step Method
2장에서 설명한 비압축성을 가정한 지배방정식은 이하와 같은 Fractional step 알고리즘에 의하여 속도 및 압력의 연성이 이루어진다. 본 연구에서는 대류항에 CIP기반의 기법이 적용되므로, 지배방정식에 유속의 공간구배에 대한 식이 필요하다. 식 (2)를 공간에 대해서 미분하면, 다음과 같은 유속의 공간구배에 대한 식을 얻게 된다.
여기서, 3차원의 경우, 아래첨자 ξ=∂()/∂ξ(ξ=x1,x2,x3), Hi는 압력항과 점성항을 나타낸다. 전체적인 계산은 아래와 같이 3단계로 이루어지게 된다.
Interface capturing step
본 계산에서처럼 자유표면이 존재하는 경우, 식 (3)을 3.3절에 설명할 THINC-WLIC기법을 사용하여 풀어, 자유수면의위치를 찾는다.
Advection step
식 (2)와 (11) 의 대류항만을 앞 절에서 설명한 USCIP법으로 풀어서 중간 속도(u*)와 구배(uξ*)를 구한다.
여기서, 첨자 n은 이전 시간스텝을 나타낸다.
Non-advection Step 1
중간속도(u*)를 이용하여, 식 (2) 중에서 압력항을 제외한 항들만을 고려한 2단계 중간속도(u**)를 구한다. 이때, 확산항은 중심차분을 이용하여 구한다.
Non-advection Step 2
속도-압력연성을 위하여, 식 (16)에 발산을 취하고 연속방정식을 적용하면 식 (17)와 같은 압력 Poisson방정식을 얻게 되며, 압력 Poisson방정식의 해법으로서는 SOR법을 채택하였다. 얻어진 압력장을 이용하여 식 (18)과 같이 다음 시간스텝(n+1)의 속도(uin+1)를 구하게 된다. 이렇게 얻어진 무발산(Divergence-free) 조건의 속도장으로부터 속도의 공간구배를 식 (19)와 같이 구한다.
3.3 THINC-WLIC기법
본 계산과 같이 다상류 유동문제를 해석하는 경우, 정확한 계면의 포착이 중요하다. 계면 포착을 위해 풀게 되는 식 (3)의 대류항에는 앞서 설명한 원형의 CIP법이나 USCIP 등을 적용할 수도 있으나, 자유표면과 같이 보다 급격한 구배의 경계를 한개의 격자의 해상도에서 좀 더 정확히 재현하기 위해서는 약간의 개선이 필요하다. 계면 포착에 사용되는 대표적인 방법인 VOF법이나 Level-set법에 사용되는 VOF 함수(혹은 밀도함수)에 Tangent 형의 변환함수를 이용하거나, 오버샷(Overshot)을 제어하기 위한 제한자(Limiter)를 사용하는 방법, 보존형 지배방정식을 사용하는 세미-라그란지안적 접근법인 CIP-CSL법, Xiao et al.(2005)에 의한 THINC법이 그 예이다. 최근에는 Yokoi(2007)가 THINC법의 확장용이성에 덧붙여 원형 이상의 정확도를 갖는 THINC-WLIC기법을 제안하였다. WLIC기법은 VOF(혹은 밀도함수)를 사용하여 계면을 포착할 때, 주어진 밀도함수 ϕij의 분포로부터 Noh와 Woodward가 제안한 SLIC(Simple line interface calculation)법(Noh and Woodward, 1976)과 유사하게 특성함수 χi,j를 구하여 자유표면을 구성하는 방법으로, SLIC과 같은 간편성을 가지면서, 면법선(Surface normal) 정보를 보다 유용하게 이용하도록 한 방법이다. WLIC기법에서는 Fig. 2와 같은 수직경계와 수평경계의 정보를 모두 사용하며, 식 (20)과 같이 면법선벡터 nij에 의해서 가중 처리된다.
위 식에서 ωx와 ωy, χx와 χy는 Fig. 2에 보이는 수직방향과 수평방향의 경계면에 대한 가중치와 특성함수이다.
본 연구에서는 식 (21)과 같은 1차원 THINC스킴을 특성함수로 채택하였으며, 차원분할(Dimensional splitting)을 통해 다차원 문제에 적용하였다.
여기서 αx와 β는 각각 경계면의 방향과 특성함수의 평활화(Smoothing)를 의미하며, αx는 아래와 같이 정의된다.
이때 1차원 법선벅터 nx,i=ϕi+1−ϕi−1이다. β는 단위 계단함수의 부드러움의 정도를 나타내며, 본 연구에서는 격자간격에서 완만한 기울기를 갖는 β=3.5를 사용하였다.
4. 수치 시뮬레이션
4.1 슬로싱 문제
개발된 계산코드의 기본 검증 문제로서 실험 및 계산결과(Yang and Kim, 2010)가 존재하는 슬로싱 문제에 대한 수치 시뮬레이션을 수행하여 결과를 비교하였다. 초기설정은 Fig. 3과 같으며, 전 방향 등간격 (Δx=0.01m)인 직각 격자가 사용되었다. 탱크의 수평 주기운동은 XH=Asin(2πt/T)을 따르며, 진폭과 주기는 각각 0.02m와 1.166s 이다. 운동효과는 격자를 움직이지 않고, 외력으로 부여하도록 하였다. 수평 및 높이 방향의 격자수는 각각 벽면에서의 속도와 압력의 경계조건으로, 각각 No-slip조건과 Neumann 조건을 사용하였다. Intel(R) Core(TM) i7-2600 (3.40 GHz) CPU와 16GB RAM을 가진 계산기를 사용하였으며, 20주기까지의 시뮬레이션에 소요된 시간은 약 510초이다.
탱크 좌측 벽면 상의 높이 0.115m 위치(P6) 및 우측 천정부분(P7)에서 측정된 압력의 시계열 데이터를 실험결과와 비교하여 Fig. 4와 5에 각각 나타내었다. Fig. 6에 나타낸 자유수면 형상을 포함하여 시뮬레이션 결과가 실험결과와 비교적 좋은 일치를 보이고 있음을 알 수 있다. 폭(y) 방향 중앙단면에서의 속도성분의 등고선(Contour)도와 경계를 표시한 Fig. 7을 보면, 기체영역에서 강한 와류를 포함한 복잡한 유동이 관찰되며, 2차원 계산에서 무시되는 폭 방향 속도성분(v) 또한 큰 값을 나타냄을 알 수 있다.
4.2 수조 내에 사각구조물이 존재하는 붕괴파 문제
지진해일과 같은 단일 대형 파도가 연안 구조물에 작용하는 충격하중을 시뮬레이션을 통해 Arnason(2005)이 수행한 실험과 비교하였다. 초기조건으로서는 Fig. 8에 보이는 것처럼 길이×폭×높이가 1.6×0.6×0.3m인 사각 수조내에 0.4×0.6×0.3m의 물기둥이 존재하며, 실험에서는 수조의 밖으로 완전배수가 불가능하여 바닥에 0.01m의 높이로 물이 채워져 있었으므로, 시뮬레이션에서도 동일하게 설정하였다. 수치 시뮬레이션에 사용된 격자는 전방향 등간격(0.01m)으로 총 576,000개의 격자가 사용되었다. 시간 간격(Δt)은 0.0005s로, 총 계산시간 3s의 시뮬레이션에 소요된 시간은 약 15,540초이다.
Fig. 9는 시뮬레이션에서 얻어진 몇몇 시간대에서의 자유표면의 형상을 나타낸다. 초기 상태에서 붕괴된 물기둥은 바닥에 위치한 물과 상호작용하며 사각구조물을 향해 전진하며, t=0.4s 일 때 구조물의 전면과 충돌하며 복잡한 형상을 나타낸다. 이 후, 수조의 오른쪽 끝 벽에 도달한 붕괴파는 반사되어 t=1.4s 부근에 기둥의 후면에 재충돌하게 된다.
이 시간에서의 구조물의 중심점을 기준으로 추출한 xz 및 yz 면에서의 압력분포, 자유표면의 위치 및 속도 벡터장을 Fig. 10에 표시하였다. 붕괴파가 구조물 및 벽에 충돌하면서 기체 영역에 강한 와류가 생성되고, 기체부분의 유동이 매우 복잡해짐을 확인할 수 있다. 공기와 물의 밀도차가 크기 때문에 그 영향은 크지 않다고 생각되나, 보다 정밀한 해석을 위해서는 이에 의한 영향을 고려해 볼 필요가 있다고 생각된다.
Fig. 10
Snapshots of air-water interface profiles with velocity vector fields and contour maps of pressure
Fig. 11은 사각구조물의 전-후면에 작용하는 파력을 Arnason(2005)의 실험과 비교한 시계열 그래프로, 실선으로 나타낸 수치 시뮬레이션 결과는 최대 충격력이 나타내는 시간에 맞추어 보정한 값이다. t=1.3~1.5s부근에서의 약간의 차이를 제외하면, 전체적으로 실험과 수치시뮬레이션 결과가 잘 일치하는 것을 알 수 있다.
5. 결 론
본 연구에서는 자유표면의 거동, 기체 및 액체의 유동과 이에 동반하여 발생하는 구조물에 대한 유체 충격 하중에 관한 문제를 해석하기 위해서 개발된 수치기법을 소개하였다. 개발된 수치기법은 다상유동의 수치해석기법의 하나인 CIP기반의 CCUP법을 근간으로 한다. 대류항과 계면 포착에 각각 사용된 USCIP법과 THINC-WLIC기법은, 기존 기법들에 비하여 수치확산이 적다는 원형 CIP법과 THINC법의 장점은 물론, 이의 확장형인 RCIP법이나 CIP-CSL계열의 다른 기법들 보다 빠른 계산시간과 적은 메모리를 사용하며, 비교적 간단한 계산으로 비선형성이 강한 이상류 문제를 정확하고 효과적으로 해석할 수 있는 장점을 가지고 있다.
개발된 수치기법을 이용하여, 사각형 탱크 내부의 슬로싱 유동과 사각구조물이 위치한 수조내에서 물기둥이 붕괴하는 경우에 대한 3차원 시뮬레이션을 수행하였다. 시뮬레이션으로부터 얻어진 자유표면의 형상과 특정위치에서의 압력의 시간변화는, 실험 결과와 비교적 좋은 일치를 보였다.
CIP기반의 기법들을 물체적합 격자나 비구조 격자에 적용하는 경우 정도가 떨어지는 것으로 알려져 있으므로, 본 계산 기법을 보다 효과적으로 사용하기 위해서는 직각 격자계를 기반으로 내/외부 물체 등을 표현할 수 있는 수치기법의 검토와 도입이 필요하다. 이를 통해, 자유수면을 포함한 복잡한 형상의 물체의 내/외부 다상유동, 나아가 물체의 운동까지 포함시키는 해석이 가능해질 것이라 기대된다.
PDF Links
PubReader
Full text via DOI
Download Citation
Print
