2.1 지배방정식

본 연구에서 제안하는 SIB법은 Lee and Mizutani(2009)와 유사하게 비압축성 점성유체를 대상으로 식(1)의 연속방정식과 식(2)의 운동방정식, 그리고 자유수면의 추적을 위해 도입하는 유체율 밀도함수 ϕ에 대한 식(3)의 이류방정식으로 구성된다.

여기서, x_i 는 데카르트 좌표에서 i = 1, 2는 x, z에 대응하며, u_i 는 i 방향의 유속성분(u, w), p는 압력, ρ는 유체의 밀도, ν는 동점성계수, Q는 소스항, t는 시간, g는 중력가속도, D_ij는 변형속도텐서, γ는 계산영역의 측면경계에 위치한 감쇠영역에서 양의 값을 갖는 감쇠계수, F_i는 SIB 법에서 고려하는 계산영역내의 구조물에 의한 i 방향의 외력항이다. 식(2)의 운동방정식은 소스항 Q를 통해 목표하는 파랑과 흐름의 발생이 가능하지만 본 연구의 목적은 가동물체형 구조물 거동해석을 위한 SIB 법의 개발이므로 Q = 0를 적용하였다.

한편, Lee and Mizutani(2009)의 IB 법의 경우 액체만을 해석대상으로 하는 단상류의 모델인 반면, 본 연구에서는 자유수면 부근에 설치되는 구조물을 정도 높게 해석하고 자유수면을 포함하는 이동경계의 효과적인 추적을 위해 Lee et al.(2011)에 의해 제안된 혼상류 모델인 단일유체모델(one-field model for immiscible two-phase fluid)을 동시에 적용하였다. 따라서, 식(2)의 유체밀도 ρ^와 동점성계수 ν^는 각 상의 유체 ϕ_i (i = 1~3; ϕ₁: 고체상태, ϕ₂: 액체상태, ϕ₃: 기체상태)의 밀도함수가 국소질량의 중심과 함께 이동하는 것으로 가정하여 다음과 같이 주어진다.

여기서, 첨자 s, w 및 a는 각각 고체, 액체, 기체의 상태를 나타내며, 유체밀도함수 ϕ_i는 ϕ₁ + ϕ₂ + ϕ₃ = 1의 관계를 갖는다. 식(1)~(3)의 지배방정식은 직교 교호격자를 적용한 유한차분법에 의해 이산방정식을 구성하였으며, 압력 연성에는 Amsden and Harlow(1970)에 의해 개발된 SMAC(simplified marker and cell) 법에 기초한 다단계(Fractional step) 법을 적용하였다.

2.2 CIP 법에 의한 수치확산의 최소화

이산화방정식은 유한한 격자점에서 구성되므로 식(6)과 같은 이류항의 계산 시 수치진동이나 확산이 발생하는 것은 잘 알려져 있다.

본 연구에서는 이류항에 대한 안정된 계산을 위해 Kim et al.(2020)의 연구에 기초하여 CIP 법을 사용하였다. Fig. 3(a)에 보인 바와 같이 각 격자점에서 초기값이 이류속도 u에 의해 ∆t시간 후 이류된 물리량의 분포형태는 Fig. 3(a)의 점선과 같은 원형 그대로 이류되어야 함에도 불구하고 Fig. 3(b)와 같이 제한된 격자점에서의 정보로 인해 이류대상 함수의 형태가 붕괴되는 것(그림 중의 파란색 실선)을 확인할 수 있다.

Fig. 3.

Principle of CIP method.

CIP 법에서는 이와 같은 물리량의 형태를 유지하면서 이류시키기 위해 식(7)과 같이 이류 물리량의 기울기(gradient)를 추가로 고려하여 이류시키는 방법으로 Fig. 3(c)와 같이 이류대상 물리량의 원형과 근사한 형태를 유지할 수 있다.

여기서, g_x, g_z는 각각 ∂f/∂x와 ∂f/∂z로 정의된다. 특히, CIP법은 인접한 격자점의 정보만으로 3차 정도의 계산이 가능하여 쌍곡선방정식의 수치해석에 있어서 수치확산이 작고 안정적인 알고리즘으로 알려져 있다(Kim et al., 2020). 인접한 격자점에 f와 g 값이 주어지면 이 격자점들 사이의 물리량 분포형태는 다음과 같은 3차 방정식으로 보간될 수 있다.

여기서, X = (x - x_i), Z = (z - z_k)이며, C_{l, m}은 다항식의 계수로 격자점 (i, k)와 인접한 격자점에서의 함수 f(x, z)와 함수의 공간기울기 g(x, z)로부터 결정할 수 있다.

2.3 구조물(고체상태)의 취급

Lee and Mizutani(2009)의 IB 법은 앞서 제시한 Fig. 2와 같이 Lagrangian point를 이용하여 유체영역을 점유하고 있는 구조물의 형상을 표현하고 구조물 경계에서 유체력을 평가하여 유체내의 구조물 영향을 Euler 격자에서 계산되는 유체장으로 피드백하는 과정을 통해 구조물과 유체의 상호작용을 모의하였다. 이와 같은 IB 법의 가장 큰 장점은 Lagrangian point를 통해 구조물의 경계를 표현하므로 임의 형상을 갖는 구조물과 유체와의 상호작용을 모의할 수 있다는 점이다. 하지만, 하지만 Lagrangian 격자와 Euler 격자간의 유체력 교환에 있어서 복잡한 보간이 필요하므로 3차원으로의 확장이 곤란하며 수치해석코드도 매우 복잡하다. 따라서, 본 연구에서는 Xiao et al.(1997)에 의해 제안된 방법인 고체상태를 강체구조물로 취급하는 방법을 적용하였다. 따라서, 물체의 운동은 이동과 회전의 조합으로 구성된다는 가정에 기초하여 해석을 수행한다.

우선 물체의 병진운동은 고체의 밀도함수 ϕ₁와 고체 밀도 ρ_s를 적분하여 계산영역 내를 점유하고 있는 물질의 총질량과 고체의 중심위치를 산정하면 다음과 같이 고체의 중심위치에서의 속도(ui¯)s를 얻을 수 있다.

여기서, (xi¯)s는 고체의 중심위치, M_s는 계산영역에 있어서 물질의 총량이다. 식(9)를 차분화하여 고체의 중심위치에서 병진속도를 구할 수 있다. 또한, 상기의 병진운동과 유사하게 고체의 회전운동에 대해서도 고체의 밀도함수 ϕ₁를 이용하여 고체의 회전축 주변의 관성모멘트 I_s를 다음과 같이 정식화 할 수 있다.

여기서, 고체에 작용하는 총 회전력(torque) Γstotal은 물체가 점유하는 각각의 격자에서 발생하는 회전력의 총합이라고 가정하면 식(11)과 같은 관계가 성립한다.

식(10) 및 식(11)의 관계를 이용하여 식(12)와 같이 각가속도를 산정할 수 있다.

따라서, 식(9)와 식(12)에 의해 계산격자 내에서의 물체의 국부이동속도 u_s를 산정할 수 있다.

2.4 SIB 법의 개요

본 연구에서 제안한 수치모델은 유체 중에 이동하는 고체입자 해석을 위해 Kashijima et al.(2001)에 의해 제안된 속도장에 대한 체적가중평균법을 Lee et al.(2011)에 의해 제안된 혼상류 모델과 결합하였다. 이하에는 본 연구에서 적용한 SIB 법을 간략히 기술한다. 기존의 IB 법은 유체입자의 운동방정식에 외력항을 도입하여 구조물의 경계면에서 산정되는 유체력으로부터 수치적으로 경계조건을 만족시키는 directingforcing 법이나 구조물 내부에 가상셀(ghost cell)을 위치시켜 보간을 통해 경계조건을 만족시키는 계산기법들이 대부분이다. 즉, IB 법은 이와 같이 구조물의 영향을 고려하기 위해 도입되는 외력항의 평가방법에 의해 구분되며 최근에는 구조물의 소성변형까지 고려할 수 있는 다양한 방법들이 제안되고 있다. 본 연구에서는 Kashijima et al.(2001)가 제안한 고체의 밀도함수로 속도장을 체적가중평균(Volume-weighted average)하는 식(13)을 적용하였다.

여기서, u_i는 흐름장에서의 속도, u_f는 국부유체속도, u_s는 국부고체속도이다.

식(13)에 의해 표현되는 흐름장의 속도 u_i는 Navier-Stokes 운동방정식을 만족하는 것으로 가정하고 계산영역 내에 존재하는 고체영역에 대해서도 고밀도의 유체로 가정하여 식(2)의 운동방정식으로부터 다음식과 같이 시간발전에 대한 가유속 u~i를 산정한다.

식(14)의 계산에서는 식(2)의 외력항인 F_i를 고려하지 않고 산정되며 계산된 가유속에 대해 F_i를 이용하여 구조물 내의 유속을 조정한다. 즉, ϕ₁ = 1이 되는 구조물셀일 경우 Fi = (us-u~i)/∆t 를 적용하고 ϕ₁ = 0이 되는 기체와 액체셀이 경우에는 F_i = 0을 적용함으로서 구조물의 이동을 수치적으로 모의할 수 있다. 또한, 0 < ϕ₁ < 1과 같이 구조물과 유체가 경계를 이루는 셀에 대해서는 Fi =ϕi (us-u~i)/∆t을 적용하면 식(14)와 유사하게 시간발전된 유체장의 속도를 식(15)와 같이 간단하게 산정할 수 있다.

식(15)에 보이는 바와 같이 본 연구에서 적용하고 있는 SIB법에 기초한 최종적인 수치결과는 FAVOR 법과 유사하게 격자내의 고체상태 점유율의 계산정도에 의존하는 것을 알 수 있다. 하지만, Lagrangian point와 Euler 격자를 혼용하는 Lee and Mizutani(2009)의 IB 법에 비해 단순한 알고리즘으로 인해 계산비용이 적고 3차원으로의 확장성이 뛰어난 장점을 지닌다.

3.1 아르키메데스 모델

본 논문에서는 새롭게 제안된 SIB 법을 수치파동수조의 적용성을 검토하기 Fig. 2에 보인바와 같이 수면상에 놓인 각주의 액체혼입깊이에 대해 수치해석결과와 아르키메데스 원리에 따른 이론 결과를 비교하였다. 아르키메데스 원리의 이론결과는 개발된 수치해석 코드의 벤치마킹으로 자주 활용되는 문제 중의 하나이다.

Fig. 4에서 보인 기체를 포함한 계산영역은 해상도 ∆x = ∆z = 1의 단위 격자를 갖는 100 × 100의 격자시스템에서 해석되었다. 각주의 높이는 H = 10∆z와 H = 15∆z를 각각 적용하였으며, 각주의 초기위치는 각주의 하단부의 위치가 정수면과 동일하게 위치하여 각주의 초기속도는 0이다. 이와 같은 계산조건 하에서 각주의 질량을 변화시켜 최대혼입깊이를 이론치와 비교한 결과를 Fig. 5에 보인다.

Fig. 4.

Configuration of Archimedes model.

Fig. 5.

Comparison between the numerical solution and the Archimedes model.

Fig. 3에서 확인되는 바와 같이 수치계산결과는 각주의 밀도에 따른 액체혼입깊이에 대한 이론치의 전체적인 경향을 성공적으로 모의하고 있음을 알 수 있다. 각주의 밀도 변화에 따른 이론값과의 오차는 수치계산에서 각주의 이동에 따라 압력과 속도의 변화로 인한 흐름장의 변화가 발생하는 반면에 이론에서는 구조물 주변의 어떠한 흐름장도 고려되지 않는 이상적인 결과에 따른 것으로 판단된다.

3.2 자유수면으로 낙하는 구조물에 대한 추적성능

본 연구의 목적은 유체 중에 운동하는 가동형 구조물의 경계면 추적을 수치적으로 모의할 수 있는 수치모델을 구축하는 것이므로 밀도가 다른 혼상의 유체(기체와 액체)에서 가동형 구조물의 추적성능을 검토하기 위해 공기중에서 자유수면으로 자유낙하하는 물체(각주)에 대한 수치모의를 진행하였다. 자유낙하 실험은 투명 아크릴로 제작된 수조에 크기 8 × 7.2 c m, 질량 m = 68 g의 물체를 수면에서 5 cm 이격 된 높이에서 자유낙하 시켰으며, 자유낙하하는 과정은 초고속 카메라(MGL 10000, TSI Inc)로 촬영하였다. 자유낙하 실험 시 각주가 공기중에 회전하지 않고 수중에 연직으로 낙하하도록 반복 실험을 수행하였다.

한편, 수치계산에서는 각 방향의 격자해상도를 ∆x = 0.4 와 ∆z = 0.2를 적용하였으며, 본 연구에서는 이동하는 물체의 경계면 추적이 목적이므로 자유수면 부근으로 계산영역을 제한하였다. Fig. 6에서부터 Fig. 8에 자유낙하하는 각주의 실험 결과와 계산된 결과를 시간경과에 따라 비교하였다. Fig. 6은 자유낙하하는 물체가 수면에 부딪히기 직전까지의 결과로 본 연구에서 제안한 SIB 법은 공기 중에서 이동하는 경계면을 잘 추적하고 있음이 확인된다.

Fig. 6.

Measured and Computed free-falling object before water impact.

Fig. 7.

Measured and Computed free-falling object at the moment of water impact.

Fig. 8.

Measured and Computed free-falling object after water impact.

Fig. 7과 Fig. 8은 자유낙하하는 각주가 수면과 접촉하여 수면아래로 진입되는 과정을 보여준다. t = 0.1103 s 이후에 각주가 수면으로 진입되면서 낙하속도가 감소하게 되고, 이러한 각주의 속도 변화는 자유표면의 파동을 형성시키는 일련의 과정을 재현하고 있음을 확인할 수 있다. 하지만, 각주의 수면진입에 따라 발생하는 자유수면의 형성패턴은 실험결과와 다소 차이가 존재한다. 이는 수치계산의 경우 단면 2차원에 대한 계산결과로 수리모형실험의 3차원적인 효과가 반영되지 못한 결과에 기인하는 것으로 판단되지만 보다 정량적인 검토가 요구된다. Fig. 8에서는 각주가 수중으로 2/3 이상이 진입하면서부터 수면에서 발생한 파동이 계산영역의 측면으로 이동하는 현상이 관찰된다. 또한, Fig. 7과 유사하게 실험에서 관찰된 각주 주변의 수면파형을 과대평가하고 있음이 확인되지만 전체적으로 각주의 수면 진입에 따라 각주의 형상과 발생되는 수면파형을 양호하게 재현하고 있다고 판단된다. 이상으로부터 본 연구에서 제안한 SIB 법은 다상의 유체가 혼합된 흐름장에서 고체의 밀도함수로 속도장을 체적가중 평균(Volume-weighted average)함으로써 유체중에 놓인 가동형 물체의 해석이 가능함을 확인하였다.