수면 아래에서 TNT(trinitrotoluene), 펜톨라이트와 같은 폭발물의 폭발로 발생한 버블은 맥동과 동시에 연직 상승하고 함정과의 상호작용에 의한 버블젯 현상 등에 의해 큰 피해를 초래한다(Zhang et al., 2011a; Zhang et al., 2011b). 따라서 함정의 구조적 안정성을 고려한 설계를 위해 버블 운동을 이해하는 것이 필요하다. 본 연구팀에서 펜톨라이트를 이용하여 수조에서 진행한 수중폭발 실험결과에 따르면(Choi et al., 2017), 폭발물이 수면에서 위치한 깊이에 따라서 폭발로 생성된 버블이 하강 혹은 상승하는 현상이 관찰되었으며, 이러한 현상은 버블과 자유수면 및 수조벽의 상호작용에 의해 결정된다고 알려져 있다(Brennen, 1995; Zhang et al., 2013). 이러한 버블과 구조물의 상호작용은 유체 내에서 회전하는 블레이드를 설계할 때 역시 필요하다. 유체의 속도 변화에 의한 압력변화로 유체 내에 공동이 발생하는데, 이러한 공동이 블레이드의 수명에 큰 영향을 미치기 때문이다(Sedlář et al., 2015; Sreedhar et al., 2017; Taskar et al., 2017; Zhang et al., 1989).

따라서 이러한 버블-구조물 상호작용에 의한 현상을 정확히 이해하기 위해서는 버블 주변의 유동을 정확히 예측하는 것이 필요하지만, 실험을 통해서 이러한 버블의 운동을 관찰하기 위해서는 ms(millisecond) 수준에서 관찰 가능한 카메라와 여러 센서 등을 포함한 가시화 장비가 필요하다는 물리적, 시간적인 한계가 존재한다(Cui et al., 2016). 이러한 한계를 극복하고자 속도 포텐셜을 이용한 이론연구와 전산 유체 시뮬레이션 연구가 진행되어 왔다. 전산 유체 시뮬레이션 방법 중 경계요소법을 이용하여 시간에 따라 변화하는 버블 거동에 대해 연구하였으며, 속도 포텐셜을 이용하여 비점성, 비회전성을 갖는 유체내의 유동장을 예측하였다(Klaseboer et al., 2005). 경계요소법을 이용하여 수치적으로 계산된 속도 포텐셜은 버블 표면과 구조물 표면에서 속도 경계조건을 만족하도록 하였다(Chahine and Perdue, 1990). 구조물과의 상호작용에 의한 효과를 연구하기 위해 구조물 부분은 Lagrangian 좌표계에서 계산을 수행하고 유체부분은 Eulerian 좌표계를 이용하였으며(Kaumuck et al., 1995; Zhang et al., 2011a), 연성 해석을 통해서 유체가 고체 구조물을 향해 이동하고 버블젯 현상을 모사할 수 있었다(Zhang et al., 2001). 이러한 컴퓨터 시뮬레이션을 이용한 해석은 정확한 계산시간에 비해 많은 계산시간을 요구한다는 단점이 있다(Zhang et al., 2011a). 버블이 항상 구 형상을 한다는 가정을 하여 속도 포텐셜을 이론적으로 계산한 기존의 연구는 적은 계산량으로 속도 포텐셜을 예측할 수 있었고 버블 표면과 강체 평판에서 경계조건을 근사적으로만 만족하며, 그 오차에 대한 체계적인 연구가 부족하다(Brennen, 1995). 속도 포텐셜을 이용한 계산은 비회전 비점성인 유동장에 대해서만 성립하지만, 해석적인 해로부터 물리적인 직관을 얻을 수 있는 장점이 존재하고, 해석적인 해를 구한 후 기존 시뮬레이션 모델의 정확도를 점검할 수 있는 테스트베드로 활용할 수 있다.

수중 폭발 혹은 공동현상으로 발생한 버블은 맥동과 병진, 2가지 독립적인 운동이 선형 중첩된 형태의 운동을 한다. 속도 포텐셜을 이용한 연구에 따르면, 무한한 유체 내에서 운동하는 버블은 두 가지 운동에 의한 속도 포텐셜을 선형적으로 더하여 전체 유동의 속도 포텐셜을 구할 수 있다. 한 방향으로 무한한 평판의 고체나 자유 표면의 경계를 갖는 유체 내에서 버블의 운동에 의한 속도 포텐셜을 계산하기 위해 기존의 연구에서는 강체(자유표면) 경계면에서의 속도 경계조건을 만족시키기 위해 강체(자유표면) 경계면 기준 대칭인 위치에 하나의 가상 버블이 있다고 가정한 후, 두 속도 포텐셜을 더하여 속도 포텐셜을 계산하였다(Brennen, 1995). 가상 버블에 의한 속도 포텐셜을 선형적으로 더하게 되면 강체(자유 표면) 경계면에서의 경계조건은 만족하지만, 가상 버블이 없을 때 만족했던 버블 표면에서의 경계조건은 가상의 속도 포텐셜이 생김에 따라 만족하지 않게 된다. 이러한 한계점을 극복하고자 이전의 연구에서는 버블과 평판사이의 거리가 버블의 반지름 보다 충분히 큰 경우에 대해 만족하는 속도 포텐셜 근사식을 Taylor 급수 전개를 통해 계산하였다(Brennen, 1995).

본 연구에서는 무한한 유체 속에서 운동하는 버블에 의한 3차원 속도 포텐셜을 이용하여 자유 표면 혹은 강체 평판과 상호작용하는 버블 주변의 속도 포텐셜을 다수의 가상 버블을 이용하여 예측한다. 병진 운동에 의한 속도 포텐셜의 경우 하나의 가상 버블을 사용한 기존 연구의 한계점을 뛰어넘기 위해 가상 버블의 개수를 늘려가며 더해가면서 강체(자유 표면) 경계면과 버블 표면 두 개의 경계조건을 만족하는 것을 보인다. 맥동에 의한 속도 포텐셜은 기존의 버블 표면에서의 오차 정도를 유지하면서 강체(자유 표면)경계면에서의 경계조건은 만족하도록 가상 버블의 개수를 늘려가며 속도 포텐셜을 계산하였다. 따라서, 본 연구에서 제시하는 속도 포텐셜을 이용한다면 버블과 강체(자유 표면)의 다양한 조건에서 해석적으로 빠르게 계산하고, 물리적인 직관을 주는 장점을 이용하여 강체와 버블 사이의 인력을 이해하는 연구가 가능하고 최종적으로 수중 폭발로 생겨난 버블이 함정에 주는 피해의 정도를 예측가능하다.

무한 유체 내에서 맥동하지 않고 병진운동만 하는 버블에 의한 속도 포텐셜은 다음과 같이 나타난다.

W는 Fig. 1에서의 x방향 속도를 나타낸다. 극 좌표계를 도입하여 r, θ를 구하면 임의의 점에서 속도 포텐셜을 계산할 수 있다. Fig. 1과 같이 버블 중심으로부터 h(=H/2)거리에 강체(자유 표면) 경계면이 존재하는 경우의 속도 포텐셜을 계산하기 위해 가상 전하 방법(Griffiths, 2013)을 사용한다. 전자기학에서 전위 경계조건을 맞추는 가상 전하 방법과 유사한 원리로, 강체(자유 표면) 경계면을 기준으로 버블의 중심에 대칭인 점에 가상 버블의 중심이 위치하게 되면 식 (2)와 같이 속도 포텐셜을 예측할 수 있으며 이는 강체(자유 표면) 경계면에서의 경계조건을 자동으로 만족하게 된다.

Φ_t는 실제 버블과 가상 버블의 속도 포텐셜의 합으로 구해진 전체 속도 포텐셜, 는 i번째 가상 버블의 병진운동에 의한 속도 포텐셜을 의미하고 r₁는 첫 번째 가상 버블의 중심으로부터 점 P까지의 거리를 나타낸다. 식 (2)의 두 번째 항 앞에 있는 부호는 경계면의 종류에 따라 결정되며 강체 평판을 경계면으로 하는 경우 +부호를, 경계면이 자유 표면인 경우는 −부호가 된다. 첫 번째 가상 버블을 도입하여 속도 포텐셜을 계산하면 강체(자유 표면) 경계면에서의 경계조건(강체 평판의 경우, 자유 표면의 경우 은 만족하지만, 가상 버블을 도입하여 생긴 로 인해 버블 표면에서의 경계조건 은 만족하지 않게 된다. 버블 표면에서의 경계조건을 만족시키기 위해서 두 번째 가상 버블을 도입한다. 버블 표면의 점 θ = 0, θ = π에서 경계조건을 만족하는 두 번째 가상 버블의 위치 x₂와 속도 𝜔₂를 이용하여 아래와 같이 식 (3)을 설정한다.

버블 표면에서 속도 경계조건을 나타내면 식 (4)와 같다.

경계조건을 적용하기 위해서 r에 대해 편미분하면 식 (5)와 같이 나타낼 수 있다.

버블 표면위의 두 점 (θ = 0, θ = π)에 적용하게 되면 이 되기 때문에 각 점에서의 r₂, r₃값을 대입하여 식 (6), 식 (7)로 식을 정리하여 나타낼 수 있다.

식 (6), 식 (7)을 x₂, 𝜔₂에 대해 연립하여 풀이하면 아래와 같이 두 번째 가상 버블의 위치와 속도를 구할 수 있다.

두 번째로 도입한 가상 버블에 의해 버블 표면에서 속도 경계조건은 정확하게 만족하지만, 강체(자유 표면) 경계면에서의 경계조건은 다시 만족하지 않게 된다. 식 (2)와 식 (8)에서 볼 수 있듯이 홀수 번째 가상 버블의 속도는 이전의 가상 버블의 속도와 부호가 반대이고 크기가 같지만, 짝수 번째 가상 버블의 속도의 크기는 직전의 홀수 번째 가상 버블에 비해 크게 감소한다. 따라서 식 (9)의 n번째 가상 버블의 위치와 거리를 이용하여 식 (9)와 같이 무한개의 가상 버블을 도입하게 되면 강체(자유 표면) 경계면과 버블 표면에서 경계조건을 모두 만족하는 속도 포텐셜을 구할 수 있다.

r_n은 n번째 가상버블의 중심으로부터 거리를 뜻한다. 기존의 연구에서는 버블과 강체(자유 표면) 경계면 사이의 거리가 버블 반지름에 비해 충분히 크다고 가정하여 한 개의 가상 버블에 의한 속도 포텐셜을 구하였는데, 이는 식 (9)에서 2번째 항(n = 1)까지 고려한 속도 포텐셜을 사용하였기 때문에 버블 표면에서 경계조건을 만족하지 않는다(Fig. 2).

본 연구에서는 가상 버블의 개수를 늘려가며 속도 포텐셜을 계산하여 각 경계조건이 어떻게 변화하는지 살펴보기로 한다. 먼저 강체 경계면을 기준으로 대칭인 버블 분포를 갖는 홀수개의 가상 버블을 더하는 경우 자동적으로 강체 경계면에서 경계조건을 만족하기 때문에 버블 표면에서의 경계조건을 살펴보면 Fig. 3과 같다. 가상 버블의 개수에 따라서 버블 표면의 한 점(r=R, θ=0)에서의 r방향 유속을 비교한 결과이며, 버블 중심과 강체 경계면사이의 거리가 버블의 반지름에 비해 클수록(h > R)더 적은 수의 가상 버블을 더하여도 경계조건을 비교적 정확하게 만족하는 속도 포텐셜을 구할 수 있다. h/R=1.2인 경우는 3개의 가상 버블을 사용할 때 대략 10% 상대오차를 갖지만, h/R=1.8의 경우는 같은 수의 가상 버블에 대해 매우 작은 오차를 나타낸다.

가상 버블의 개수를 늘려가면서 버블 표면에서 속도 경계조건을 살펴보면, 약 9개의 가상 버블을 이용하여 계산할 경우 버블과 강체 경계면이 매우 가까운 h/R=1.2인 경우에도 표면에서 경계조건을 만족한다. 따라서 9개의 가상 버블을 이용하여 점(r=R, θ=0)뿐만 아니라 임의의 θ에 대해서 버블 표면에서 반지름 방향 속도를 그려보면 Fig. 4에서 볼 수 있듯이 표면 모든 점에서 경계조건을 만족하는 것을 볼 수 있다.

짝수개의 가상 버블을 도입하여 속도 포텐셜을 계산할 경우 버블 표면에서 자동적으로 경계조건을 만족하고 강체 경계면에서의 경우 그 오차가 가상 버블의 개수에 따라 달라진다. 강체 경계면의 한 점(x=h, y=0)에서 v_x를 계산하면 Fig. 5와 같다. 버블의 중심이 강체 경계면으로부터 충분히 먼 경우 비교적 적은 수의 가상 버블을 사용하여 속도를 예측하여도 그 오차가 상대적으로 작으며, 8~10개의 가상 버블을 도입하여 계산할 경우, 강체 경계면으로부터 버블이 충분히 가까워도 무시할만한 오차를 보이며 경계조건을 만족한다.

강체 경계면에 수직한 방향의 속도 값을 강체 경계면의 위치에 대해 계산하면 Fig. 6과 같다. h/R=1.2인 경우에 대해 10개의 가상 버블(N^img=10)을 더할 경우, 그 오차는 매우 작아진다. 앞서 가상 버블의 개수가 홀수 일 때의 결과와 종합하여 보면, 가상 버블의 개수가 10개 이상이면 그 숫자가 홀, 짝 종류에 상관없이 두 개의 경계조건을 모두 만족하는 것을 알 수 있고 만약 무한개의 가상 버블을 더하게 되면 그 오차는 0으로 수렴하게 된다. 10개의 가상 버블을 이용하여 h=1.2R인 경우 경계조건을 모두 만족하는 속도 포텐셜을 xy평면에서 그려보면 Fig. 7로 나타낸다. 속도 포텐셜을 공간에 대해 수치 미분하여 화살표로 유동장을 가시화 하였고 모든 변수는 W,R로 나누어 무차원화 하였다.

본 연구에서 제시하는 방법의 오차를 고찰하기 위해서 주어진 계산 자원(Intel Core i7-6700 CPU @ 3.40GHz)을 이용하여 병진운동하는 버블 주변 3차원 속도 포텐셜을 계산하고 그 결과의 오차를 정리하면 Table 1과 같다. 벽면에서의 이론값은 0이기 때문에 벽면에서의 속도를 W로 나눈 값을 그대로 나타냈고, 버블 표면에서의 오차는 이론값에 대해 상대 오차를 나타냈다. 홀수, 짝수에 상관없이 10개 이상의 가상 버블을 이용하여 본 연구에서 제시하는 3차원 속도 포텐셜을 계산할 경우 짧은 계산 시간으로 높은 정확도를 갖는 것을 알 수 있다.

버블의 맥동운동에 의한 속도 포텐셜은 식 (10)과 같다.

은 버블 반지름의 시간에 대한 변화율을 의미한다. 강체(자유 표면) 경계면(x=h)에서 경계조건을 만족시키기 위해 강체(자유 표면) 경계면을 기준으로 대칭인 위치에 첫 번째 가상 버블을 도입하면 속도 포텐셜은 식 (11)로 나타난다.

두 번째 항 앞에 있는 부호는 강체 경계면의 경우 −, 자유 표면 경계인 경우 +로 계산된다. 맥동에 의한 속도 포텐셜의 경우, 병진에 의한 속도 포텐셜과 달리 버블표면에서 경계조건을 정확하게 만족시키는 2번째 가상 버블의 위치와 반지름의 시간에 대한 변화율을 구하는 것이 수학적으로 어렵기 때문에 경계조건을 오차 (R/H)² 정도까지 만족하는 속도 포텐셜을 계산한다. 2번째 가상 버블에 의한 속도 포텐셜 ()을 더하고 버블 표면에서 경계조건을 적용하면 식 (12)와 같다.

코사인 제2법칙을 이용하여 r₂를 계산한 뒤, 식 (12)에 대입하여 정리하면 아래와 같다.

이전의 연구에서 식 (14)의 우변을 Taylor 근사하고 (R/H)³이상의 고차항은 값이 작다고 가정하여 식을 계산하였고, 그 결과, 버블 표면에서 주어진 오차까지 경계조건을 만족하는 2번째 가상 버블의 맥동에 의한 속도 포텐셜을 식 (15)로 나타낼 수 있다(Brennen, 1995).

2번째 가상 버블의 중심이 1번째 버블과 일치하기 때문에 조건 r₂=r, θ₂=θ을 사용하여 식 (15)에 나타냈다. 이전의 연구에서는 식 (15)에 표현된 2번째 가상 버블까지 고려하여 속도 포텐셜을 계산하였는데, 2번째 가상 버블의 도입으로 강체 경계면에서 경계조건을 만족하지 않는다. 본 연구에서는 이러한 기존연구의 한계를 해결하고자 추가적인 가상 버블을 도입하여 계산하기로 한다. 먼저 강체 경계면을 기준으로 대칭인 위치에 3번째 가상 버블을 두면, 강체 경계면에서의 경계조건을 만족하고 버블 표면에서 만족하지 않는다. 버블 표면에서 기존의 오차정도 (R/H)²를 갖으며 동시에 강체 경계면에서의 경계조건을 만족시키기 위해서 앞서 설명한 방법을 이용하여 4번째 가상 버블의 위치와 버블 반지름의 시간에 대한 변화율을 계산할 수 있다.

식 (16)의 가상 버블에 의한 속도 포텐셜은 병진운동에 의한 속도 포텐셜과 식의 형태가 유사하기 때문에 맥동에 의한 속도 포텐셜을 구할 때 사용했던 방법을 이용하여 아래와 같이 무한급수의 형태로 표현이 가능하다.

맥동운동에 의한 속도 포텐셜은 함수 형태가 달라 병진운동에 의한 속도 포텐셜과 달리 두 경계조건을 모두 정확히 만족하는 함수를 계산하기 어렵다. 그 결과, Fig. 8에 나타나있듯이, 무한개의 가상 버블을 중첩하여 속도 포텐셜을 계산해도 버블 중심과 강체 경계면의 거리에 따라 버블 표면에서 경계조건을 만족하지 않는다.

버블 표면에서 경계조건 오차는 (R/H)³정도이기 때문에 h/R 값이 커질수록 오차의 크기가 작아지기 때문에 버블 표면에서 경계조건을 더 정확하게 만족하게 된다. h/R>1.8인 경우 버블 표면에서 10% 미만의 상대 오차를 갖는 반면, h/R=1.2인 경우 충분히 많은 가상 버블을 중첩하여 계산하여도 30% 이상의 상대오차를 갖는다. h=2R일 때 9개의 가상 버블을 중첩하여 표면에서의 속도를 Fig. 9와 같이 나타낼 수 있다. 버블 반지름의 변화율로 무차원화된 속도 변수는 버블 표면에서 대략 5% 미만의 오차를 갖으며, θ=0에서 최대 오차를 갖는다.

버블의 중심과 강체 경계면의 거리에 상관없이 많은 가상 버블을 도입하면 강체 경계면에서 경계조건을 만족하는 속도 포텐셜을 구할 수 있기 때문에 작은 h/R 값에 대해서도 경계 조건을 만족한다. 대략 11개의 가상 버블을 사용하게 되면, 강체 경계면에서는 h/R=1.2인 경우에도 경계조건을 만족한다. 강체 경계면(x=h)에서 경계조건을 가상 버블의 개수(N^img)를 증가시키면서 계산해보면 Fig. 10처럼 오차가 감소하는 것을 알 수 있다. 11개의 가상 버블을 중첩시켜서 계산할 경우 강체 경계면에서 10^-5 수준의 오차를 갖는다(Fig. 11).

11개의 가상 버블을 중첩시켜서 x=h=2R의 위치에 강체 평판이 있는 경우에 대해 계산한 무차원 속도 포텐셜과 유동장을 가시화하면 Fig. 12와 같다. 강체 평판에서 x방향으로의 속도가 0으로 경계조건을 만족하며, 속도 포텐셜을 갖는 이상 유체이기 때문에 강체(자유 표면) 경계면에서의 y(x)방향 속도는 유한한 값을 갖는다.

맥동 운동하는 버블에 의한 속도 포텐셜 계산 결과를 정리하면 Table 2와 같다. 계산 소요 시간은 병진 운동의 경우와 마찬가지로 대략 0.5초 정도가 소요되며, 만족해야하는 경계조건에 대한 오차도 작게 나타났다. 벽면에서의 경계조건은 비교적 정확하게 만족하는 반면, 버블 표면에서의 오차는 상대적으로 크게 나타났다. 이러한 오차는 앞서 설명한 Taylor 급수 전개를 통해 얻은 계산의 오차이며 그 값은 10⁻² 정도로 나타났다.

본 논문에서는 가상 버블을 이용하여 맥동, 병진운동하며 강체(자유 표면)경계면과 상호작용하는 버블 주변 유동의 속도 포텐셜을 이론적으로 예측하였다. 이전의 연구에서는 2개(병진) 혹은 3개(맥동)의 가상 버블을 사용하여 속도 포텐셜을 계산하였기 때문에 병진운동에 의한 속도 포텐셜의 경우 버블의 표면에서, 맥동운동에 의한 속도 포텐셜의 경우 강체(자유 표면) 경계면에서 경계조건을 만족하지 않았다. 병진 운동에 의한 속도 포텐셜은 가상 버블의 개수를 대략 10개까지 늘려가면서 포텐셜을 계산하면 버블 표면과 강체(자유 표면) 경계면에서 경계조건을 동시에 만족하는 속도 포텐셜을 계산 할 수 있었다. 맥동운동의 경우는 속도 포텐셜 함수 형태가 병진 운동에 의한 속도 포텐셜과 달라서 직접적인 가상 버블을 도입하기는 어려움이 있었다. 이에 본 연구에서는 이전의 연구에서 제시한 방법을 이용하여 버블 표면에서 (R/H)³ 정도의 오차를 갖는 속도 포텐셜을 사용하였고, 더 나아가 기존의 연구에서 강체(자유 표면) 경계면에서 경계조건을 만족하지 않는 한계점을 가상 버블을 추가적으로 도입하여 극복하였다. 후속연구로써, 맥동하는 버블의 경계조건을 더 정확하게 만족하는 해석적인 풀이 방법을 개발하고, 그렇게 구해진 속도 포텐셜을 이용하여 강체(자유 표면) 경계면에 의해 버블에 가해지는 인력과 척력을 계산하는 연구를 계획 중에 있다. 본 연구에서 제시한 속도 포텐셜은 버블과 고체의 상호작용으로 발생하는 인력을 이용하여 버블에 의한 함정의 피해 가능성을 예측하는 연구에 적용 가능할 것이라 예상되며 또한, 수조에서 수중 폭발 실험을 하는 경우 자유 표면과 강체(바닥)으로부터 인력과 척력을 받기 때문에 발생하는 버블의 상승 혹은 하강을 예측하는 연구가 가능하다.