계산비용의 지속적인 감소와 비약적인 계산성능의 향상에 기인하여 해안 및 해양공학 분야에서도 유체의 운동방정식을 공간 및 시간에 대한 이산화와 다양한 수치해석적인 기법을 도입하여 공학적인 문제를 해결하려는 전산유체역학(CFD)적인 시도가 활발하게 이루어지고 있다. 초기의 전산유체역학은 단상유동(single phase flow)을 중심으로 발전해 왔지만, 최근 들어서는 자유수면을 포함하는 파동장과 같이 단상의 경계가 시간발전에 따라 지속적으로 변화하는 경우나 액상과 기상이 혼합되는 문제에 있어서는 다상유동(multi phase flow) 문제를 적용하는 사례가 급격하게 증가하고 있다. 특히, 파동장과 같은 자유수면의 문제를 취급하는데 있어서는 혼합되지 않는 액상과 기상의 비압축성 뉴턴유체를 고려한 혼상류 모델이 적용되는 경우가 많다.

일반적으로 서로 혼합되지 않는 유체(multiple immiscible fluids)는 개별 유체가 갖는 속성의 급격한 변화(또는 불연속성)가 존재하므로 각각의 영역으로 구분될 수 있다. 즉, 유체간의 경계면(자유수면)은 하나의 유체(기체) 밀도가 다른 유체(액체)와 비교하여 무시할 정도로 작다는 극단적인 형태의 경계면으로 취급할 수 있다. 일반적으로 혼상류 모델은 이상에서 언급한 각상의 경계면에 대한 시간기반 거동추적이 필수적이며, 궁극적으로는 유체간 경계면(자유수면)의 추적성능이 계산의 정도를 좌우할 수 있다.

일반적으로 경계면 추적기법에는 유한요소법, 경계적합 좌표계법, 라그랑주(Lagrange) 법과 같이 계산격자를 경계면 변형에 추종하면서 분석하는 방법과 계산격자의 재구성을 필요로 하지 않는 방법인 Marker and Cell(MAC) 법(Welch et al., 1966), Level-Set 법(Osher and Sethian, 1988), Volume of fluid(VOF) 법(Hirt and Nichols, 1981) 등으로 구분할 수 있다. 이와 같은 다양한 경계면 추적기법 중에서 오일러 격자체계를 이용하여 각각의 계산격자에 대한 유체의 점유율에 주목한 VOF 법이 가장 널리 활용되어 왔으며, 최근 들어서는 Level-set 법을 적용한 계산 예도 지속적으로 증가하고 있다. 대표적인 공중사용허가서 라이선스 기반의 CFD 프로그램인 OpenFOAM(Jacobsen et al., 2012; Higuera et al., 2013)에서도 자유수면에 대한 추적기법으로 VOF 법을 적용하고 있다. VOF-type의 경계면 추적법은 유속(流束, flux)을 도출함에 있어서 반드시 자유수면에 대한 형상정보가 필요하며, 이러한 자유수면에 대한 형상정보는 계산범위 및 격자해상도에 비례하여 계산부하를 증가시킨다. 한편, Lee et al.(2012) 은 자유수면에 대한 형상정보가 요구되지 않는 대수적 경계 추적기법인 THINC 법(Xiao et al., 2005)을 파동장의 해석에 적용하여 계산정도를 평가하였다. Keshavarzi et al.(2013)는 2차원에서 부력에 의해 상승하는 기포의 형태와 변형을 VOF 법과 Level-set 법을 결합한 Coupled Level Set and VOF (CLSVOF) 법을 이용하여 해석하고, 기존의 VOF 법과의 해석성능을 비교하였다. 그 결과, 기존의 VOF 법은 기포의 운동을 정확하게 모사하기 위해서는 CLSVOF 법에 비해 보다 상세한 격자구성이 요구되며, 반대로 CLSVOF 법은 경계면이 분리 및 재결합되는 곳에서는 정확하게 경계면을 추적할 수 없음을 지적하였다. Jabbari et al.(2014)은 비뉴턴 유체의 혼상류를 대상으로 VOF 법을 적용하여 서로 다른 세 가지 경계면 재구축 방법(Geometric Reconstruction Scheme(GRS), High Resolution Interface Capturing(HRIC), Compressive Interface Capturing Scheme for Arbitrary Meshes(CICSAM))의 적용성을 검토하여 CICSAM을 이용한 VOF 법이 낮은 계산비용으로 비뉴턴 유체의 경계면 추적에 가장 적합함을 보였다. Xie et al.(2017)은 복잡한 경계를 갖는 혼상류해석을 위해 비구조 격자체계에서 유한체적법에 기초한 THINC 법을 적용한 UMTHINC(unstructured multi-dimensional THINC) 를 제안하고 기존의 VOF법에 비해 높은 경계면 추적성능을 나타냄을 증명하였다. 이외에도 다양한 분야에서 경계면 해석을 위한 연구들이 다수 수행되어 왔으나 주로 경계면 추적 기법의 성능향상을 위한 연구이며, 다양한 경계면 추적기법의 경계면 추적성능을 정량적으로 상호 비교한 연구예는 매우 제한적이다.

본 연구에서는 상용 및 오픈소스 기반의 다양한 CFD 수치해석코드에 적용되어 활용되고 있는 대표적인 오일러 격자 체계에 기반한 VOF-type의 경계면 추적기법들의 이류성능을 평가하기 위해 Zalesak(1979)에 의해 제안되어 이류계산의 벤치마킹에 자주 적용되는 Zalesak disk를 대상으로 하였다. 특히, 기존의 전통적인 VOF-type의 경계면 추적법 및 이류계산에서 발생하는 수치확산을 최소화 하기 위해 수치유속(numerical flux)을 제어하는 FCT 법(Boris and Book, 1973)의 효용성을 평가하고, 더불어 CIP 법(Yabe and Aoki, 1991)을 활용한 자유수면 추적성능의 가능성을 고찰하였다.

일반적인 VOF-type 경계면 추적기법은 해석하고자 하는 대상유체(해안공학분야에서는 액체에 해당됨)에 의해 점유되는 셀에서는 1의 값을, 그렇지 않은 셀에서는 0의 값을 갖는 단위체적당 유체의 점유체적비를 나타내는 다음의 유체율 함수(이하 VOF 함수)를 적용한다.

여기서, x는 위치벡터이다.

식(1)의 VOF 함수가 0과 1의 사이의 분수를 갖는 셀에서 경계면의 정보를 제공하므로 VOF 함수를 지시함수(indicator function)라고 불린다. 즉, VOF-type 경계자유수면 추적기법은 각각의 계산격자에 하나의 숫자로 표현되는 최소화된 저장정보로 경계면 추적이 가능함에 따라 Hirt and Nichols(1981)의 SOLA-VOF 법이 제안된 이후에 다양한 공학분야에서 요구되는 경계면 추적해석에 적용되어 왔다. 잘 알려진 바와 같이 VOF 함수 F의 시간의존은 다음의 이류방정식에 의해 지배된다.

여기서, t는 시간, u는 대상 유체의 속도벡터이며 VOF 함수가 대상 유체와 함께 이동함을 의미한다.

따라서, 식(2)의 이류방정식을 유한차분근사(Finite-difference approximation)하여 시간발전에 따른 경계면의 추적이 가능하다. 하지만, 일반적인 차분근사를 적용하는 경우에는 저차정도 근사에 따른 수치확산 및 고차정도 근사에 따른 수치진동에 의해 경계면 정의를 손실할 수 있다. 이에 Hirt and Nichols(1981)는 VOF 함수의 시간발전을 보다 정도 높게 계산하기 위하여 VOF 함수가 컬러함수(color function)임을 활용한 donor-acceptor method(Johnson, 1970)를 적용하였다. Hirt and Nichols(1981)가 적용한 donor-acceptor method에 의한 이류면에서의 VOF 함수 F_AD의 기본적인 개념을 Fig. 1에 보인다.

Fig. 1에서 아래첨자 D와 A는 각각 doner 셀과 accpetor 셀을 의미하며, u_i+1/2 은 이류면에서의 이류속도, △x는 셀의 크기, △t는 시간간격, V = u_i+1/2△t로 이류속도에 의해 이동한 거리를 각각 나타낸다. Fig. 1에 보인 바와 같이 이류면에서의 VOF 함수 F_AD는 doner 셀과 acceptor 셀의 조건에 따라 결정되며, 이류플럭스 δF는 다음과 같이 계산된다.

여기서, Min은 donor 셀이 보유 유체량보다 더 많은 유체가 이류되는 것을 방지하고, Max는 acceptor 셀의 보유 기체량보다 더 많은 기체가 이류되는 것을 방지한다(Kim et al., 2002).

본 연구에서는 이상의 Hirt and Nichols(1981)의 SOLA-VOF 법에 의한 경계면의 추적정도를 검토하기 위해 Fig. 2와 같이 Zalesak(1979)에 의해 제안된 예제(이하에서는 Zalesak disk라고 칭함)에 적용하여 이류성능을 검토했다.

Zalesak disk 초기 프로파일은 식(4)와 같이 반경 0.17m와 폭 0.05 m, 길이 0.27m의 슬롯을 갖고 있다.

여기서, R = {(x-0.5)²+(x-0.75)²}^1/2 이며, 계산영역의 격자 간격은 △x = △y = 0.01 m, 격자해상도는 100 × 100이며, 계산 시간 간격 △t = 1.0 sec으로 설정하였다. 또한, 이류속도는 다음의 식(5)를 적용하였다.

여기서, x_c와 y_c는 각각 계산영역의 중심좌표이며 α는 CFL을 조정하기 위한 계수로 본 연구에서는 CFL_MAX 가 각각 0.314, 0.157, 0.079가 되도록 설정하였다.

Figs. 3~5는 CFL 조건에 따른 SOLA-VOF 법에 의한 Zalesak disk의 경계면 추적성능을 나타내며 Figs. 3~5(a)는 초기조건의 Zalesak disk를 반시계방향으로 1/2 회전 시킨 경우를, Figs. 3~5(b)는 초기조건의 위치까지 1회전 시킨 경우를 각각 나타낸다. 계산결과로부터 확인되는 바와 같이 SOLA-VOF 법에 의한 경계면의 추적성능은 초기의 Zalesak disk의 형상을 정도 높게 재현하지 못함을 확인할 수 있다. 특히, VOF 함수의 이류계산 시 이류속도에 의한 완전한 플럭스 수송이 발생하지 못함에 따라 Zalesak disk의 이동궤적에 따라 계산영역 내에 작은 잔상(droplet)들이 발견된다. 이러한 계산영역 내에 발생하는 잔상들은 CFL 조건이 가장 높은 Fig. 3에서 가장 두드러지고 CFL 조건이 감소함에 따라 다소 개선된 계산결과를 보이지만 여전히 계산영역에는 잔상이 존재함을 알 수 있다. CFL = 0.079의 경우는 CFD 계산 시 계산비용이 급격하게 증가하는 결과를 가져오므로 실제로 적용할 수 없음에도 불구하고 Fig. 5에서 보이는 바와 같이 잔상을 제거할 수 없다. 더욱이, 작은 잔상들은 경계면의 추적 성능 저하와 더불어 질량보존을 담보할 수 없다. 이상의 결과에서 확인되는 바와 같이 Hirt and Nichols(1981)가 제안한 본래의 SOLA-VOF 법은 계단함수인 VOF 함수의 특성을 유지하기 위해 적용한 donor-acceptor 스킴만으로는 정도 높은 경계면 추적을 기대할 수 없음을 의미한다.

Hirt and Nichols(1981)가 제안한 SOLA-VOF 법 이후로 Torrey et al.(1985, 1987)은 표면장력과 벽면점착을 추가적으로 고려한 NASA-VOF code를 제안하였다. 또한, 셀에서 유체가 차지하는 면적비를 이용하여 계산격자와 일치하지 않는 벽면경계를 고려하는 fractional volume of fluid method를 코드에 추가로 삽입하였다.

한편, 일반적인 이류계산에 있어서 오차를 제어하기 위하여 Maximum and minimum Bounds preserving(MmB) 스킴(Wu and Yang, 1989)이나 limiter function을 일반적으로 적용한다(Lee et al., 2012). 이러한 VOF 함수의 이류관점에서 NASA-VOF 코드의 가장 큰 특징은 VOF 함수의 값을 초과하는 값을 갖는 표면셀에 대하여 표면셀의 체적, 기체셀의 체적 및 액체셀의 체적을 이용하여 표면셀의 값을 조정한다는 점이다. 실제 코드에 적용된 표면셀의 VOF 함수 조정에 관한 계산코드를 아래에 제시한다.

if(f[i][k]>1.0) {

f[iec][kec]=f[iec][kec] + VS * (f[i][k]-1.0) / VE;

f[i][k]=1.0;

}

if(f[i][k]<0.0) {

f[ifc][kfc]=f[ifc][kfc] + VS * f[i][k] / VF;

f[i][k]=0.0;

}

여기서, VS는 표면셀의 체적, VE는 기체셀의 체적, VF는 액체셀의 체적을 각각 나타내며, [iec][kec]는 기체셀의 셀번호, [ifc][kfc]는 액체셀의 셀번호를 의미한다.

Fig. 6은 Fig. 4와 동일한 CFL 조건을 적용한 경우에 대한 NASA-VOF 법의 Zalesak disk 형상재현정도를 나타낸다. 그림으로부터 확인되는 바와 같이 SOLA-VOF 법을 적용한 Fig. 4에 비해 NASA-VOF 법을 적용한 Fig. 6의 결과가 Zalesak disk의 잔상이 현저하게 감소하였음을 확인할 수 있다. 즉, 유체점유율 함수에 대한 limiter function을 적용한 SOLA-VOF 법에 비해 유체율 함수가 이상값을 갖는 표면셀을 대상으로 주변셀의 체적을 이용하여 표면셀의 값을 조정하는 NASA-VOF 법이 유체율 함수에 대한 이류계산을 보다 정도높게 계산할 수 있음을 확인하였다.

FCT(Flux-corrected-transport) 법은 오일러 유한차분 격자체계에서 일반적인 연속방정식을 수치적으로 안정하게 해석하기 위해 Boris and Book(1973)에 의해 제안된 이후로 다양한 이류계산에 FCT 법이 활용 또는 개선되어 왔다(Zalesak, 1979; Murray, 1997). FCT 법은 언더슈트(undershoot) 또는 오버슈트(overshoot)가 발생할 가능성이 높은 영역(이류대상의 물리량의 구배가 강한영역)에서 유속(flux)을 제한하도록 설계된 알고리즘으로 일반적으로 2단계로 나뉜다. 즉, 이류계산 단계와 확산방지 또는 수정단계의 두 가지 주요단계로 구성된다.

본 연구에서는 다음과 같은 순서로 FCT 법을 수행하였다.

1) 저차정도의 스킴과 고차정도의 스킴을 이용한 셀 경계에서의 유속 Fi+1/2L 및 Fi+1/2H를 계산

2) 수치확산 방지플럭스를 A_i+1/2 = Fi+1/2H - Fi+1/2L와 같이 정의

3) 저차정도에 기초한 유속을 이용하여 이류계산을 수행, Witd=Win-∆xt-1Fi+1/2H-Fi+1/2L

4) 수치확산방지 유속의 한계치(임계치) 계산, AiC=Ci+1/2 Ai+1/2, 0≤CCi+1/2≤1

5) 수치확산방지 유속을 적용한 이류방정식 계산, Win+1=Witd-∆xt-1Ai+1/2C-Ai-1/2C

Fig. 7은 Fig. 6 및 Fig. 4와 동일한 CFL 조건을 적용한 경우에 대한 FCT-VOF 법의 Zalesak disk 형상재현정도를 나타낸다. Fig. 3의 SOLA-VOF 법과 Fig. 6의 NASA-VOF 법을 적용한 VOF 함수의 이류계산에서 나타났던 잔상이 FCT-VOF 법에서는 수치확산을 최소화 하기 위해 수치유속을 제한함으로써 제거되었음을 확인할 수 있다.

CIP(Constrained Interpolation Profile or Cubic Interpolated Pseudo-particle) 법은 쌍곡선형 방정식을 해결하기 위해 Takewaki et al.(1985)에 의해 제안된 방법으로, 그 후로 다양한 형태의 CIP 법이 개발되어 왔다(Takewaki and Yabe, 1987; Yabe and Aoki, 1991; Yabe et al., 2001, 2004). CIP 법의 기본적인 개념은 이류계산에 있어서 격자점의 정보(함수)와 함께 이류대상 함수의 기울기 정보를 동시에 계산하는 방법이다. 식(6)의 1차원 이류방정식을 고려하면, CIP 법의 이류계산에는 식(7)과 같이 함수의 기울기를 고려한 방정식을 동시에 고려한다.

여기서, G는 함수 F의 공간미분 G = ∂F/∂x를 나타낸다. 2개의 격자점에 F와 G 값이 주어지면 이 격자점들 사이의 프로파일은 식(8)과 같은 3차 방정식으로 보간될 수 있다.

여기서, a_i, b_i, c_i 및 d_i는 격자점 i에서의 3차 방정식의 계수이다. 한편, 2차원의 경우는 식(8)의 3차 방정식은 식(9)와 같이 확장된다.

여기서, X =(x - x_i), Y =(y - y_i)이다. A_type CIP 법의 경우는 식(9)로부터 다음의 식(10)과 같이 10개의 계수를 갖는 3차 방정식으로 된다.

식(10)의 계수 C_l,m은 격자점 (i, j)와 인접한 격자점인 (i+1, j), (i, j + 1) 및 (i+1, j + 1)에서의 함수 F(x, y), 그리고 격자점 (i, j), (i+1, j), (i, j + 1) 및 (i+1, j + 1)에서의 함수의 각방향으로의 기울기 ∂F(x, y)/∂x, ∂F(x, y)/∂y를 적용하여 결정할 수 있다.

Fig. 8과 Fig. 9는 Fig. 3 및 Fig. 4와 동일한 CFL = 0.314, 0.157를 각각 적용한 경우에 대한 A_type CIP 법에 의한 Zalesak disk 형상 재현정도를 나타낸다.

Fig. 8로부터 확인되는 바와 같이 SOLA-VOF 법이나 NASA-VOF 법에서 나타났던 작은 잔상들은 관찰되지 않지만 수치확산에 의해 Zalesak disk의 경계면 형상이 완만하게 재현되며, disk 주변에서는 언더슈트가 발생한다. CFL 조건이 낮은 Fig. 9의 경우도 Fig. 8과 같이 유사한 재현결과를 보이지만 오히려 full revolution 후에는 Fig. 8에 비해 Zalesak disk의 형상 재현정도가 감소함을 확인할 수 있다.

CIP 법은 기존의 VOF-type의 경계면 추적기법과는 상이하게 VOF 함수의 이류계산 이후에 경계면에 대한 형상의 구축절차가 필요하지 않으며 계산 복잡성에서 상당히 경제적이다. 하지만, Fig. 8과 Fig. 9에서와 같이 급격한 불연속성이 나타나는 부분에서 수치확산 혹은 수치진동이 발생할 수 있다. 이를 해결하기 위해 VOF 함수의 값이 대부분 0과 1 부근에 집중되어 있는 점을 고려하여 VOF 함수의 변환을 통해 국지적인 공간해상도를 향상시킬 수 있다. 따라서, VOF 함수 F를 대신하여 변환함수 H(F)를 적용하면 식(2)의 이류방정식은 다음과 같이 변환될 수 있다.

여기서, 변환함수 H(F)는 tangent 함수를 이용하여 다음과 같이 정의된다.

여기서, ε은 작은 양의 정수이다. 또한, 변환함수 H(F)의 이류계산 이후에 VOF 함수는 식(12)의 역변환인 식(13)으로부터 구할 수 있다.

Fig. 10으로부터 VOF 함수의 tangent 변환을 통해 국지적인 공간해상도를 향상시킨 결과를 확인할 수 있다. 즉, 원래의 VOF 함수 F는 0에서 1까지 급격하게 변화하는 반면에 변환함수 H(F)는 0과 1 근처에 집중되어 있음을 확인할 수 있다.

Fig. 11은 VOF 함수의 tangent 변환을 통한 CIP 법의 계산결과를 나타낸다. Fig. 11과 동일한 계산조건인 Fig. 8의 결과와 비교하면 tangent 변환을 적용한 Fig. 11 의 경우가 Zalesak disk 주변의 공간해상도를 국지적으로 향상시킴에 따라 재현정도가 향상되었음을 확인할 수 있다. 또한, SOLA-VOF 법과 NASA-VOF 법에서 발생하는 작은 잔상들이 발생하지 않음을 알 수 있다.

본 연구에서 적용한 대표적인 VOF-type 경계추적기법인 SOLA-VOF, NASA-VOF, FCT-VOF 및 CIP 법에 의한 오차율을 비교하였다. 오차의 계산은 식(14)를 적용하였다.

여기서, Fi,jO는 셀(i, j)에서 VOF 함수의 초기조건을 Fi,jT는 이류계산을 수행한 후의 VOF 함수를 나타낸다.

본 연구에서 적용한 경계추적기법에 따른 오차를 Table 1에 보인다. SOLA-VOF 법, NASA-VOF 법 및 FCT-VOF 법과 같이 이류계산 후에 경계면에 대한 형상의 재구축이 필요한 Geometrical-type의 경계추적법의 경우는 CFL 조건이 감소함에 따라 계산오차도 감소한다. 반면에 이류계산 후에 경계면에 대한 형상의 재구축이 불필요한 CIP 계산의 경우에는 CFL 조건이 감소함에도 불구하고 오히려 계산오차가 증가하는 경우가 발생할 수 있으므로, CIP 법을 선택하는 경우에 있어서는 격자해상도를 신중하게 선택할 필요가 있음을 시사한다. 또한, CIP 법의 경우에는 VOF 함수의 tangent 변환을 통해 국지적인 공간해상도를 향상시켜 이류계산을 수행함으로써 오차를 대폭 감소시킬 수 있음을 확인하였다.

최근 들어, 계산기의 성능이 비약적으로 향상됨에 따라 CFD를 이용하여 다양한 공학문제를 해결하려는 시도가 증가하고 있으며, 해안공학분야에서도 쇄파대에서의 메커니즘, 난류구조의 해석 및 지형변동에 이르기까지 적용범위가 확대되고 있다. 특히, 공중사용허가서 라이선스 기반의 다양한 CFD Tool이 개발되고 배포됨에 따라 사용자는 모델에 삽입된 계산 스킴을 선택하여야 하는 경우가 종종 발생하고 있다. 이러한 CFD 모델별로 제공되는 다양한 계산 스킴의 옵션 중에서 해석대상의 목적에 따른 적절한 경계면 추적기법의 선택은 최종적인 계산결과에 직접적인 영향을 미칠 수 있으므로 매우 신중해야 할 필요가 있다. 현재까지 배포된 오일러 격자기반의 CFD Tool은 경계면 추적기법을 포함하고 있으며, 그 중에서도 VOF-Type의 경계면 추적기법이 주를 이루고 있다.

본 연구에서는 이러한 VOF-Type의 경계면 추적기법의 성능을 Zalesak disk를 통해 비교·분석하였다. 그 결과를 요약하면 다음과 같다.

1) SOLA-VOF 법은 VOF 함수가 컬러함수임을 고려하여 donor-acceptor 스킴을 적용하여 이류계산을 수행하고 있지만, donor-acceptor 스킴만으로는 정도 높은 경계면 추적을 기대하기 어렵다.

2) NASA-VOF 법은 이류계산시 발생하는 이상값(outlier)을 갖는 표면셀에 대해 인접한 액체셀과 기체셀의 체적비를 통해 표면셀의 값을 조정함에 따라 SOLA-VOF 법에 비해 경계면의 추적성능이 향상됨을 확인하였다.

3) FCT 법은 고차정도와 저차정도를 동시에 이용하여 수치확산의 방지를 도모하는 기법으로 본 연구에서 논의한 VOF-type의 경계면 추적기법 중에서 가장 높은 추적성능을 보였다.

4) CIP 법은 VOF 함수의 이류계산 이후에 경계면에 대한 형상의 구축절차가 필요하지 않으며 계산 복잡성에서 상당히 경제적이지만 경계면 주변에서 수치진동과 수치확산이 발생하였다. 이러한 수치진동과 수치확산은 tangent 변환을 통해 국지적인 공간해상도를 향상시킨 결과 경계면의 추적성능이 증가함을 확인할 수 있었다.

본 연구에서 이류계산의 대상으로 적용한 Zalesak disk는 경계면이 극단적으로 변화하는 경우로 실제 파동장과 반드시 대응될 수 없다. 하지만, 다양한 수치해석코드에 적용되는 자유수면의 추적기법을 선택함에 있어서 중요한 기초자료로 활용될 것으로 기대된다.