전산유체역학을 이용한 기포장막의 수리학적 현상 재현 연구

1. 서 론

수중 압축공기의 방출로 이루어진 기포장막(Bubble curtain)의 주변 유동장에 관한 연구는 한 세기 동안 많은 연구자들의 관심분야로 여겨져 왔으며 이러한 메커니즘을 활용한 다양한 응용분야가 존재한다.

대표적으로 해수면 흐름 생성으로 항만 또는 항구지역을 내습하는 파랑으로부터 보호하기 위한 연구(Taylor[1955]; Bulson[1961]; Bulson[1963]; Bulson[1968]; Brevik[1976]; Jang et al.[2020]), 음용수(Drinking water) 저수지에 조류(Algae) 성장을 방지하기 위한 미세기포의 주입(Schladow[1992]), 유조선 사고 또는 해상공사로 인한 유류 오염 발생 시 유막이 퍼지는 것을 막거나 기름으로 인한 피해로부터 해안 서식지 보호(Fannelop and Sjoen[1980]), 조선소의 갑문수로(Lock-channel)를 통한 염수 침입 저감(Abraham and Burgh [1962]), 해양 및 하천 빙결 방지 그리고 물위를 떠다니는 부유체를 이동시키는 방안(Riess and Fannelop[1995]) 등이 있다.

또한 기포장막의 공학적 이론에 관한 초기 연구로서 Taylor[1955]는 수중 방출되는 기포의 난류 흐름과 열원에서 연기의 거동 이론과의 유사점을 연계하였고 추후 이 이론은 현상학적 분석으로 연구가 진행되었다(Brevik and Kristiansen[2002]). Kobus[1968]는 기포장막의 거동에 대해 상세 분석 모델을 적용하여 수심 2.0, 4.3, 10.4 m에서 실험을 수행하고 기포장막의 수직 속도 프로파일(Profile)을 가우스 분포(Gaussian distribution)로 추정하였다. 그리고 Bulson[1961]; 1963]; [1968]은 기포장막의 외부 흐름의 대표적인 대규모 실험 연구로서 최대 수심 10 m에서 기포 공급원에 의하여 유도된 수평 흐름의 크기를 측정하고 그 속도와 생성 두께에 대한 경험식을 유도하였다.

하지만 기포장막의 합리적이고 경제적인 작동을 위한 전제 조건은 상승하는 기포에 의하여 유도되는 기포장막 주변의 유동과 기포 공급원의 설치 매개변수 사이의 정량적 관계를 정립하는 것이 요구된다. 따라서 많은 연구자들이 기포장막의 생성과 주변 유동장에 미치는 영향을 전산유체역학(Computational fluid dynamical method)을 통해 관련 연구를 수행하였으나 현재까지 기포장막 외부 난류 영역에서의 재현 정확도는 다소 미흡한 실정이다. 더구나 수중 방출되는 압축공기의 수치해석 시뮬레이션을 정도 높게 재현하기 위해서는 수치해석의 경제적인 해석을 위한 유효 격자 구성, 기체와 액체 각 상(Phase)의 운동, 기포에 대한 난류의 분산 효과 그리고 각 상의 양방향 운동량 및 에너지의 전달이 검토되어야 한다.

따라서 본 연구에서는 전산유체역학 상용프로그램인 FLOW-3D를 이용하여 Fig. 1과 같이 수중 압축공기의 방출로 유도되는 기포장막의 수리학적 현상을 재현하기 위한 수치해석을 수행하고자 하였다. 특히 실제 수리현상의 재현성 향상을 위해 수치해석 모델에서 요구되는 경험계수를 모델 민감도 분석을 통해 도출하고 선행 연구자의 기존 실험 결과와 비교함으로써 해석 모델의 재현성을 검증하고자 하였다. 최종적으로 구축된 수치해석 모델을 이용하여 기포 방출량에 따른 기포장막의 거동 특성을 분석하고 기포장막 주변 유동장에 대한 수리학적 특성을 기포 방출량과 최대 유속의 관계식으로 산출함으로써 기포장막으로 유도된 외부 흐름장의 실제 현장 적용을 위한 공학적인 기초 설계인자를 제안하고자 한다.

Fig. 1.

The bubble curtain geometry and flow variables.

2. 재료 및 방법

2.1 선행연구

Kobus[1968]는 기포장막에 대한 상세 분석 모델을 적용하여 기포 공급원(Air source)의 깊이 2.0, 4.3 m 조건에서 모형실험을 수행하였다. 수직 속도 프로파일을 가우스 분포로 분석함으로써 수중 방출된 기포로 이루어진 기포장막의 너비는 가상의 기포 공급원으로부터 수직 높이(z+z₀)에 따라 선형적으로 증가함을 주장하였다. 여기서, z는 기포 공급원으로부터 수직거리, z₀는 실제 기포 공급원 아래 가상 공급원(Virtual source)의 깊이이다(Fig. 1).

이후 Ditmars and Cederwall[1974]은 단상 부력 기포 기둥 이론을 바탕으로 등온 상태에서 기포의 팽창과 슬립 속도를 함께 고려하였다. 수직 유체 속도 및 기포로 인한 물과 기포의 혼합물의 밀도 손실에서 Kobus의 실험결과와 동일한 프로파일을 가정함으로써 기포 상승으로 인한 연행을 평균 중심선의 속도에 비례한다고 가정하여 연행률(Rate of entrainment)을 다음과 같이 정의하였다.

$\[ \frac{d{Q}_w}{dz}=2\alpha w_m \]$

(1)

여기서, Q_w는 물의 유량이고 연행계수(Entrainment parameter) α는 0.1의 값으로 일정하다고 가정하였다.

위 식 (1)은 운동방정식과 함께 수치 적분된 두 개의 적분 방정식으로부터 산출된다. 그리고 깊이 감소에 따른 기포의 팽창은 기포장막의 부력을 증가시키고 축 대칭 기포장막의 경우 내부 코어(Core)에만 기포가 존재하는 것으로 나타난다. 기포장막을 형성하는 중요한 매개변수로서 기포 코어 너비와 기포장막의 전체 폭 사이의 비율(Schmidt number, λ로 표기)에 대해 Fannelrop and Sjeren[1980]과 Milgram[1983]은 축 대칭 기포장막의 경우 λ=0.6~0.8의 범위를 가지며 실험 결과에서 연행률 α가 가스 체적 유량에 의존하고 공급원으로부터 거리에 따라 달라질 수 있음을 보여주었다. 기포장막의 기포 팽창은 다방성 법칙을 따르는 것으로 가정하여 다음과 같이 기포의 유입으로 인한 물의 밀도 변화로 정의할 수 있다.

$\[ {\mathrm{\rho }}_a\left(z\right)={\mathrm{\rho }}_a\left(0\right){\left[\frac{p\left(z\right)}{p\left(0\right)}\right]}^{1/n} \]$

(2)

여기서, ρ_a는 기포의 밀도, p(z)는 압력 깊이이다. 일정한 수온의 경우, 상승하는 기포의 팽창은 등온조건만을 고려하며 이때 다방지수(Polytropic index) n은 1의 값을 가진다.

Fannelop et al.[1991]이 제안한 수직 유동은 수평 흐름의 척도(Scaling) 매개변수와 표면 흐름의 영향 범위를 계산하는데 사용할 수 있는데 수직 속도의 시간 평균과 기포장막 중심의 밀도 분포를 다음과 같이 가우스 분포로 나타낸다.

$\[ w\left(x,z\right)=w_m\left(z\right)\exp \left(\frac{-x^2}{b^2}\right) \]$

(3)

$\[ {\mathrm{\rho }}_w-\mathrm{\rho }\left(x,z\right)=\left[{\mathrm{\rho }}_w-\mathrm{\rho }\left(z\right)\right]\exp \left(\frac{-x^2}{{\lambda }^2{b}^2}\right) \]$

(4)

여기서, w(x, z)는 x, z 위치에서의 상승 속도, w_m(z)은 공급원으로부터 높이 z에서의 기포장막의 중심선 속도, b는 기포장막의 폭, ρ_w는 물의 밀도, ρ(x, z)는 x, z 위치에서의 밀도, ρ(z)는 공급원으로부터 높이 z에서의 밀도이다.

자유 표면에서 기포장막 내부의 상승 흐름의 편향으로 수평 흐름이 유도된다. 이 경우 흐름이 수면을 따라 기포장막 양옆으로 전달될 때 표면에서의 최소한의 에너지 손실을 가정한다(Andreopoulos et al.[1986]). 이를 통해 수중에서 방출되는 기포의 공급량과 표면에서 수평 흐름 사이의 관계에 관한 연구는 Taylor[1955]; Bulson[1961]; Abraham and Burgh[1962]; Kurihara[1965]에 의해서 실험적 및 이론적 연구가 수행되었다. Taylor[1955]가 제시한 삼각형 분포를 가정하여 수평 흐름 패턴은 Fig. 1과 같으며 이 프로파일은 실제로 계측된 속도 프로파일에 가깝다. 수중에서의 기포 방출량에 따른 수평 흐름에 대한 기존 연구자들의 실험적 관계식은 Table 1에 정리하였다.

Table 1.

The formula for horizontal flow induced by the bubble curtain

2.2 수치해석 모델

본 연구에서 수중 압축공기의 방출로 유도되는 기포장막의 수리학적 현상을 재현하기 위한 수치해석 모델로 상용 코드 FLOW-3D 모델을 사용하였다.

사용된 모델의 RANS(Reynolds Averaged Navier-Stokes) 방정식은 RNG k−ε 난류모델과 결합된 유한 체적 방법(Finite Volume method)을 사용하여 수치적으로 계산되며 난류 운동에너지 k와 난류 소산율 ε은 각각의 수송방정식(Transport equation)에 의하여 모델링된다. RANS 방정식과 연속방정식을 기반으로 난류 영역을 해결하고 급격한 변동 영향을 최소화하였다. 또한 수치해석 모델에서는 각 유체의 경계 및 해석 영역의 자유 표면(Free surface)을 재현하기 위하여 VOF(Volume of fluid) 방법을, 수중에 연행된 공기를 고려하고 계산하기 위하여 난류 공기 연행(Turbulent air entrainment) 및 Drift-flux 모델을 적용한다. 여기서, 공기 연행 모델은 계산 격자보다 작은 기포의 연행 예측을, Drift-flux 모델은 지배방정식에 항력 성분을 추가하여 반송된 위상(Carrier phase)에 분산된(Dispersed) 공기 영향을 고려한다. 이는 계산 격자의 증가로 인해 물리현상(상세 계산 격자를 요구하는 소규모의 유체 역학 문제)의 재현에서 수치해석 모델을 적용하기 어려운 경우에 사용함으로써 경제적인 시뮬레이션이 가능케 한다.

2.2.1 지배방정식

본 연구에 적용된 수치해석 모델의 기본방정식은 질량 및 운동량 보존을 반영하는 유체역학에 기초하며 비압축성 뉴턴 유체를 고려한 질량 및 RANS 방정식은 다음과 같다.

$\[ \mathrm{\rho }\frac{\partial u_i}{\partial x_i}=0 \]$

(5)

$\[ \mathrm{\rho }\frac{\partial u_i}{\partial t}+\mathrm{\rho }{u}_j\frac{\partial u_j}{\partial x_j}=-\frac{\partial }{\partial x_i}\left(-p{\delta }_{ij}\right)+\mathrm{\mu }\frac{\partial }{\partial x_j}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)+g \]$

(6)

여기서, ρ는 유체의 밀도, u_i는 직교좌표계에서 각 방향(i = 1,2,3; x, y, z)으로 유속 벡터, t는 시간, g는 중력가속도, p는 압력, μ는 점성계수(dynamic viscosity) 그리고 δ_ij는 Kronecker delta로 i = j 조건에서 1, i ≠ j 조건에서 0이다.

본 연구의 CFD 해석에 사용된 RANS 방정식은 유동의 시간 평균적 성질(평균 속도, 평균 압력, 평균 응력 등)을 계산한다. 하지만 RANS 방정식으로 난류를 계산하기 위해서는 레이놀즈 응력을 예측해야 한다. 시간 평균과정에서 발생한 변수(레이놀즈 응력 성분)를 포함한 시스템을 수학적으로 닫기(closed) 위해서는 추가적인 방정식이 필요하기 때문이다. 본 연구에서는 RNG k−ε 난류 모델을 적용하여 레이놀즈 응력을 계산한다.

2.2.2 난류 모델

난류란 안정화하려는 점성력이 불충분할 때 발생하는 불안정하고 혼란한 운동이다. 이는 레이놀즈 수가 클 때 유동 내에서 발생하는 자연적 불안정성은 줄어들지 않고 여러 크기의 와류의 생성으로 나타난다. 그러나 컴퓨터 메모리와 해석 시간 제약으로 인해 일반적으로 난류 변동의 모든 현상을 명확히 수치해석하는 것은 가능하지 않다. 그러므로 평균 유동 특성에서의 난류 효과를 기술하는 단순화된 모델을 선택할 수 밖에 없다.

본 연구에서는 수중 공급원에서 압축공기의 유입으로 발생하는 난류를 계산하기 위하여 RNG k−ε 모델을 적용하였는데, 이는 Navier-stokes 방정식에 재정규화(Renormalization) 기법을 적용하여 더 작은 규모의 운동 효과를 나타낼 수 있다(Yakhot et al.[1992]).

본 연구에서 난류해석에 적용된 RNG k−ε 모델은 통계적 방법으로 난류 운동에너지(k)나 소산율(ε)과 같은 난류량들에 대하여 평균화된 방정식의 유도로 계산된다. 차원해석을 통해 운동방정식을 난류수송방정식과 연결하는데 필요한 난류 동점성 계수(v_T)는 다음과 같이 정의된다.

$\[ {\nu }_T=C_{\mathrm{\mu }}\frac{k^2}{\mathrm{\epsilon }} \]$

(7)

여기서, C_μ는 0.085이다. 난류 모델의 특정 수치문제는 ε 값을 낮게 제한하는 것이다. 운동방정식의 계산과정에서 0에 가까운 ε을 산출할 경우 난류 운동에너지도 0에 가까워야 한다. 수치상의 이유로 그렇지 않을 경우 비물리적으로 큰 난류 동점성 계수가 발생된다. 이러한 잠재적인 수치문제를 해결하기 위하여 최대 난류 확산 거리에서 최소 허용 난류 소산을 다음과 같이 정의하여 처리한다.

$\[ {\mathrm{\epsilon }}_{T,\mathrm{m}\mathrm{i}\mathrm{n}}=C_{\mathrm{\mu }}\sqrt{\frac{3}{2}}\frac{k^{3/2}}{TLEN} \]$

(8)

여기서, TLEN은 최대 난류 확산 거리이다. 이것은 경험적 결과로 특정값을 지정할 수도 있고 수치해석 프로그램의 계산(computed) 조건을 이용하여 시간과 공간이 다른 TLEN을 계산할 수 있다. 이러한 계산 과정에서 난류 길이 및 시간 규모의 하한은 Kolmogorov 규모를 기반으로 하고 상한은 급속왜곡이론(Rapid distortion theory)을 기반으로 한다(Isfahani and Brethour[2009]). 난류 길이 및 시간 규모의 상한과 하한은 다음과 같이 정의한다.

$\[ \begin{array}{cc}{L}_{T,\mathrm{m}\mathrm{i}\mathrm{n}}=70{\nu }^{3/4}{\mathrm{\epsilon }}^{-1/4},& L_{T,\mathrm{m}\mathrm{a}\mathrm{x}}=\frac{0.86}{C_{\mathrm{\mu }}}\frac{\sqrt{k}}{S}\end{array} \]$

(9)

$\[ \begin{array}{cc}{L}_{T,\mathrm{m}\mathrm{i}\mathrm{n}}=6\sqrt{\frac{\nu }{\mathrm{\epsilon }}},& L_{T,\mathrm{m}\mathrm{a}\mathrm{x}}=\frac{0.35}{C_{\mathrm{\mu }}}\frac{1}{S}\end{array} \]$

(10)

여기서, v는 동점성 계수이고 S는 레이놀즈 응력의 평균 변형률의 크기이다. 위 식에서 주어진 한계에 따른 난류 길이는 식 (8)에서 TLEN을 대체한다.

TLEN은 난류 점성 영역의 길이 범위를 제어하는 역할을 하는데 본 연구에서는 기포장막의 흔들림 현상이 없고 중심축에 대칭적인 상승 기포를 가정하여 속도 및 밀도 프로파일의 재현이 선행연구 결과와 가장 일치하도록 각 계산 시간별 최대 난류 길이를 시간과 공간의 함수로 산출하는 계산(computed) 조건으로 설정하였다.

2.2.3 기포연행 및 Drift-flux 모델

수중에서의 기포-물 상호작용인 공기 연행은 유체요소와 관련된 에너지 밀도(P_d)가 단위 체적당 난류 운동에너지(P_t)보다 작을 때 발생한다. 이때 P_t는 흐름을 불안정하게 만드는 섭동 성분을 나타내며 단위 시간당 연행된 공기의 부피(dV)는 다음과 같이 정의된다.

$\[ dV=C_{air}{A}_s\sqrt{2\frac{P_t-P_d}{\mathrm{\rho }}} \]$

(11)

여기서, C_air는 보정 매개변수이고 각 격자의 표면적이다. 이때 공기량 dV는 혼합물의 거시적 밀도에 영향을 미치고 흐름의 팽창을 생성하는 것을 고려하여 혼합물의 밀도는 다음과 같이 계산된다.

$\[ \mathrm{\rho }=\left(1-C\right){\mathrm{\rho }}_w+C{\mathrm{\rho }}_a \]$

(12)

여기서, C는 모든 격자에서 계산된 연행된 공기 비율이다. 또한 공기가 연행되어 특정 직경의 기포를 형성한다고 가정한다. 따라서 기포는 물의 움직임에 저항하는 항력을 생성하며 이때 Drift-flux 모델을 적용하여 물과 공기 상(phase) 사이의 상대 속도가 계산된다. 유체-입자(particle), 유체-기포(air bubble) 그리고 밀도가 다른 유체-유체의 혼합물로 구성된 유동장에서 각 상의 성분들은 각기 다른 유동 속도를 가지는데 이러한 속도 차이는 밀도 차이로 인해 불균일한 체적력을 가지기 때문에 발생한다.

Drift-flux 모델에서는 각 성분 사이의 상대 운동은 입자와 같은 분산요소보다는 연속체로 근사될 수 있다. 즉, 분산요소의 운동 및 상호작용의 추적을 계산할 필요가 없으므로 계산효율을 높일 수 있다. 드리프트 근사에서 상대 속도는 유동이 하나는 연속상이고 다른 하나는 분산상이며, 이는 불연속적이고 연속인 상에 의해 둘러싸여 있다고 가정하고 두 성분 유체의 비압축성 유동에 대해 ∇⋅u = 0를 정의한다. 혼합물을 구성하는 두 성분의 체적률은 F_c과 F_d로 표시하며, 연속상에서의 속도와 분산상의 속도를 각각 u_c, u_d로 정의한다. 여기서, 아래첨자 c는 연속상(continuous phase), d는 분산상(dispersed phase)의 성분을 나타낸다.

상대 속도는 각 상의 속도에 종속하므로 이때 항력은 부유상의 체적율에 대한 정보를 가져야 한다. 부유상이 같은 크기인 입자들로 구성되어있고 단위 체적당 n개가 있다고 추정하면, 이때의 상대속도는 다음과 같이 정의될 수 있다(Brethour and Hirt[2009]).

$\[ u_r=\left(\frac{V_p}{K_p}\right)\frac{F\left({\mathrm{\rho }}_c-{\mathrm{\rho }}_d\right)}{\overline{\mathrm{\rho }}}\nabla P \]$

(13)

여기서, V_p = (1 – F)/n는 단일 입자의 체적, K_p는 연속유체 내의 속력 |u_r|을 가지는 단일 입자에 대한 항력계수, $$ \overline{\mathrm{\rho }}=F{\mathrm{\rho }}_c+\left(1-F\right){\mathrm{\rho }}_d$$는 체적가중 밀도이다. U가 연속 유체 내에서 움직이는 입자의 상대속도의 크기라면 이때 단일 입자에 대한 항력계수는 다음과 같다.

$\[ K_p=\frac{1}{2}{A}_p{\mathrm{\rho }}_c\left(C_{D}U+\frac{12{\mathrm{\mu }}_c}{{\mathrm{\rho }}_c{R}_p}\right) \]$

(14)

여기서, C_D는 기포의 항력계수, R_p는 평균 입자 반경이며, A_p는 구체 형상으로 가정되는 입자의 단면적이다.

분산된 물질들의 체적율이 충분히 작지 않을 경우 성분간의 모멘텀 교환 계산을 위한 단일 입자 항력의 사용은 적합하지 않다. 따라서 입자간의 상호작용을 고려하기 위해 Richardson-Zaki 상관관계를 적용한다. 이 상관식은 드리프트 속도를 다음 계산식으로 나타낼 수 있다.

$\[ u_r^{eff}=u_r\cdot \mathrm{m}\mathrm{a}\mathrm{x}{\left(0.5,f\right)}^{\zeta } \]$

(15)

여기서, 지수 ζ(=C_rz ζ₀)는 Richardson-Zaki 계수(C_rz)와 ζ₀의 곱으로 정의된다. 조정된 유효값은 모든 Drift-flux를 계산하는데 적용한다.

수중에서 기포의 움직임은 공기와 물의 밀도 차이로 인한 부력에 의하며 공기의 밀도가 물의 밀도에 비하여 무시할 수 있는 수준이면 공기의 부력은 압력변화량 ∂P/∂z의 함수로 정의된다. 일반적으로 부력은 기포의 직경에 따른 부피 차이로 달라지며 상승속도에 따라 특정 저항을 생성한다. 그리고 기포의 변형은 구형을 복원하는 힘과 표면장력에 의하여 상쇄되는 외부 힘에 기인한다.

수중에서의 기포 운동은 다음과 같은 무차원수를 기반으로 정적인 힘의 균형을 설명할 수 있다. 웨버수(Weber number)는 표면장력과 유체 관성력의 상대적인 힘의 크기를 나타내며 다음과 같이 정의된다.

$\[ We=\frac{{\mathrm{\rho }}_{c}U{R}_p}{\sigma } \]$

(16)

여기서, σ는 기포 경계면에서의 표면장력으로 0.073 kg/sec²의 값을 가진다. 이와 유사한 이론이 Levich[1962]에 의하여 제안되었으며 입자의 내부 압력은 변형된 입자의 모세관 압력(Capillary pressure)과 균형을 이루며 다음과 같이 난류응력으로 정의된다.

$\[ We=\frac{\mathrm{\tau }{R}_p}{\sigma }{\left(\frac{{\mathrm{\rho }}_d}{{\mathrm{\rho }}_c}\right)}^{1/3} \]$

(17)

위 식은 균일한 등방성 난류를 입자 크기의 척도에 대한 난류 응력($$ \left(\mathrm{\tau }=\overline{{\mathrm{\rho }}_c{u}^{\text{'}2}}\right)$$)과 난류속도 성분(u')으로 표기할 수 있다(Batchelor[1951]).

그리고 카필러리 수(Capillary number)는 서로 섞이지 않는 두 유체의 점도와 표면장력의 상대적인 힘의 크기를 나타내는 무차원수이며 다음과 같이 정의된다.

$\[ Ca=\frac{{\mathrm{\mu }}_{c}U}{\sigma } \]$

(18)

여기서, μ_c는 연속상(물)의 동점성 계수로 0.001 kg/m/sec의 값을 가진다. 본 연구에서는 수중에서 속도를 가지는 기포의 변형과 관련하여 임계 웨버수(We_crit) 및 카필러리 수(Ca_crit)를 매개변수로 수치해석 재현성 평가를 수행하였다.

2.2.4 계산 영역 및 격자 생성

수치계산 격자의 생성처리는 CFD 모델링에서 많은 주의가 요구되는 주요 단계이다. 일반적으로 격자의 수와 크기는 수치해석 시뮬레이션의 중요한 기준이며 정확도를 제한한다. 결과 및 계산 시간, 임계 격자 크기를 결정하는 적절한 방법은 상대적으로 큰 격자로 시작하여 원하는 효율에 도달할 때까지 줄여나가는 것이다. 그 이후 추가적인 격자의 크기 감소는 결과에 영향을 주지 않는다.

본 연구에서의 계산 영역은 Fig. 2와 같이 기포장막의 유체 흐름을 분석하기 위하여 2차원 격자로 구성하였다. 기포 공급원에 가까울수록 격자 크기를 감소시켰으며 이것이 결과의 정확도에 영향을 미치지 않도록 계산에 소요되는 시간을 고려하였다. 격자 크기에 따른 기포장막 중심선 속도의 비교는 Fig. 3에 나타냈는데 그 결과, 격자의 크기가 매우 커지지 않는다면 중심선 속도에서의 영향은 크지 않은 것으로 나타났다. 본 연구에서는 기포장막의 중심에서의 속도프로파일 분포 계산을 고려하여 기포 공급원 주변에서의 Δx, Δz 방향 격자크기를 0.01 m로 최종 선정하였다. 결과적으로 본 연구에서 사용된 총 격자수는 88,220개, 격자 크기는 0.010<Δx<0.120 m, Δy = Δz = 0.010 m이다. Δx 격자크기의 경우 효율적인 수치해석을 수행하기 위해 가변격자를 적용하였으며 기포 공급원 기준으로 ±0.5 m 영역의 격자를 더욱 조밀하게 구성하였다.

Fig. 2.

Information on mesh geometry and boundary conditions at this study.

Fig. 3.

Centerline velocity (wm) versus distance above air source (z) for different grid size at Qa = 180 LPM.

경계 조건은 모든 실험에 동일하게 적용하였는데 X_min 및 X_max 경계는 유출 조건(outflow), Y_min와 Y_max에는 대칭 조건(Symmetric), Z_min은 벽면 조건(Wall)으로 Non-slip 조건을, Z_max는 압력 조건(Pressure, 1 atm)으로 설정하였다.

3. 결과 및 고찰

3.1 모델 경험계수 및 실험안 선정

본 연구에서는 수치해석 모델의 민감도 분석을 위해 Kobus[1968]의 수리 실험 결과와 비교 검증하였다. Kobus[1968]의 실험은 폭 1.0 m, 길이 10.0 m, 수심 2.0 m의 수조에 홀 크기 0.1 cm, 홀 간격 10.0 cm의 다공관을 설치하여 기포 공급량 180, 372, 600 LPM(liter per minute)의 변화에 따라 기포장막 중심에서 유속의 시간 평균 및 공급원으로부터 특정 높이(0.8, 1.1, 1.4 m)에서의 유속 프로파일을 제공하며 수직 방향 유속분포 프로파일의 폭(b)은 가상공급원으로부터 거리가 멀어질수록 선형적으로 증가하는 것으로 나타났다. 실험 결과의 민감도 해석은 기포장막의 중심선 속도(w_m) 크기와 중심으로부터 x 방향의 수직 유속 분포 프로파일을 비교하였다.

각 경험계수의 변화에 따른 민감도 해석을 위해 Table 2와 같이 실험안을 구성하여 수치해석 시뮬레이션을 수행하였다. 기포 연행 및 Drift-flux 모델을 적용한 수치해석에서 고려해야 하는 경험계수는 각 격자에 기포의 연행율을 결정하는 연행계수(C_air), 구체 형상으로 가정하는 기포의 항력계수(C_D), 초기 기포 크기(d₁), 드리프트 흐름에서 상대 속도 크기에 기여하는 Richardson-zaki 상관방정식의 상관 계수(C_r2), 임계 웨버수(We_crit) 및 임계 카필러리 수(Ca_crit)이다. 실험안을 살펴보면 S01~S28의 경우 공기배출량 180 LPM이 수중에서 방출될 때의 각 모델 매개변수의 민감도 영향을, S29~S34은 공기배출량이 372 LPM 조건일 때 임계 카필러리 수의 영향을, S35~S40은 공기배출량 600 LPM 조건일 때 임계 웨버수의 영향을 평가하였다.

Table 2.

Numerical model variables for sensitivity analysis

3.2 모델 경험계수의 민감도 분석 결과

민감도 분석 결과에 따르면 연행계수(C_air)의 값이 0.1~0.8 범위, 초기 기포 크기(d₁)가 0.005~0.020 m 범위에서 그 계수값을 변화시키면서 계산하였지만 기포장막의 중심선 속도 및 수직 유속 프로파일에서 그 변화 영향이 나타나지 않았다. Fig. 4와 Fig. 5는 기포방출량 180 LPM 조건에서 연행계수와 초기 기포 크기변화에 따른 중심선 속도를 나타내며 연행계수 및 초기 기포 크기의 변화는 기포장막의 생성에 영향을 미치지 않는 것으로 나타났다.

Fig. 4.

Centerline velocity (wm) versus distance above air source (z) for different air entrainment (Cair) at Qa = 180 LPM.

Fig. 5.

Centerline velocity (wm) versus distance above air source (z) for different initial drop dia. (d1) at Qa = 180 LPM.

또한 기포의 항력계수(C_D)가 0.05~2.00 범위를 가지는 경우 중심선 속도의 크기에 영향을 미쳤으나 유속 프로파일에는 영향이 나타나지 않았다. Fig. 6에 항력계수 변화에 따른 중심선 속도를 나타내었으며 이때 적용된 경험계수는 Table 2의 S04~S09의 조건을 적용하였다. Fig. 6을 살펴보면 C_D 값이 증가할수록 중심선 속도가 커지는 경향이 나타났다. 이는 식 (14)에서 단일 입자의 항력계수 K_p에 영향을 미치고 항력계수의 값이 커질수록 상대속도의 크기가 작아지면서 나타나는 현상으로 생각된다. 그리고 Fig. 7에 항력계수 변화에 따른 공급원으로부터 높이 1.40 m에서의 수직 유속 프로파일을 제시하였는데 프로파일의 항력계수의 변화는 프로파일 형상에 영향을 미치지 않는 것으로 나타났다.

Fig. 6.

Centerline velocity (wm) versus distance above air source (z) for different drag coefficient (CD) at Qa = 180 LPM.

Fig. 7.

Vertical velocity profiles versus distance from centerline (x) for different Cd at z = 1.40 m, Qa = 180 LPM.

다음으로 입자간의 상호작용을 계산하기 위하여 도입된 Richardson-Zaki 상관방정식의 상관 계수(C_rz)는 0.01~10.0의 범위로 중심선 속도 및 프로파일의 영향 정도를 살펴보았다. Fig. 8에 C_rz 변화에 따른 중심선 속도를 나타내었으며 이때 적용된 경험계수는 Table 2의 S12~S16의 조건을 적용하였다. C_rz의 값이 작아질수록 중심선 속도 크기가 공급원과 가까운 지점에서 작아지는 경향이 나타났으며 유속 프로파일에는 영향을 미치지 않았다.

Fig. 8.

Centerline velocity (wm) versus distance above air source (z) for different Richardson-Zaki coefficient (Crz) at Qa = 180 LPM.

외력이 작용할 때 기포의 변형과 관련된 무차원수 중 웨버수(We_crit) 및 카필러리 수(Ca_crit)의 임계값 변화에 따른 중심선 속도 및 유속 프로파일의 변화를 각각 Fig. 9와 Fig. 10에 나타내었다. Fig. 9에서 적용된 각 경험계수 조건은 Table 2의 S29~S34, Fig. 10에 적용된 경험계수 조건은 S35∼S40을 적용하였다. We_crit 및 Ca_crit의 변화를 살펴보면 기포 방출량(Q_a)이 커질수록 각 계수의 임계값이 커지는 것을 알 수 있다. We_crit 및 Ca_crit값이 커질수록 중심선 속도가 낮아졌으며 이는 각각의 임계값이 클수록 낮은 상대속도에서 기포 입자의 분해를 유발하여 중심선 속도의 감소로 나타나는 것으로 판단된다. 특히 Ca_crit의 경우 일정 임계값보다 클 때에는 중심선 속도의 변화가 나타나지 않는다. 그리고 특정 Ca_crit 값에서 We_crit는 Fig. 11과 같이 중심선 속도의 최대값을 결정하는 매개변수로서 나타낼 수 있으며 이를 이용하여 각각의 기포방출량 조건에서 목표 중심선 속도를 계산하는 We_crit 값을 산출할 수 있다. 그리고 We_crit 및 Ca_crit 값의 변화는 유속 프로파일에는 영향을 미치지 않았다.

Fig. 9.

Centerline velocity (wm) versus distance above air source (z) for different critical capillary number (Cacrit) at Qa = 372 LPM.

Fig. 10.

Centerline velocity (wm) versus distance above air source (z) for different critical weber number (Wecrit) at Qa = 600 LPM.

Fig. 11.

Maximum Centerline velocity (wm,max) for different critical weber number (Wecrit) at Qa = 372 LPM.

3.3 RNG k−ε 점성 난류 모델 경험계수

최대 난류 확산 거리인 TLEN의 경우 Table 2의 S21~S26의 조건으로 0.002~0.200 m의 범위 및 각 계산 시간별 최대 난류 길이를 시간과 공간의 함수로 산출하는 조건(computed 조건)을 함께 고려하여 중심선 속도 및 유속 프로파일의 영향을 살펴보았는데 각 경험계수는 S21~S26의 조건을 적용하였다. 그 결과는 Fig. 12과 Fig. 13에 제시하였는데 중심선 속도 및 유속 프로파일에서 영향이 모두 나타나는 것을 확인할 수 있다. TLEN 값이 작을수록 점성 영향으로 인해 상승 기포의 흔들림(Wandering) 현상이 나타났으며 TLEN 값이 공급원 깊이의 7%보다 크거나 computed 조건을 적용하였을 때 기포장막의 흔들림이 발생하지 않고 Kobus[1968]의 실험 결과와 유사하게 유속 프로파일이 재현되었다.

Fig. 12.

Centerline velocity (wm) versus distance above air source (z) for different maximum turbulent mixing length (TLEN) at Qa = 180 LPM.

Fig. 13.

Vertical velocity profiles versus distance from centerline (x) for different TLEN at z = 1.40 m, Qa = 180 LPM.

3.4 최종 경험계수 및 모델 검증

수중에서 압축기포의 방출로 생성된 기포장막의 수치해석 시뮬레이션 재현을 위해 여러 모델 경험계수들의 민감도 분석을 수행하고 모델에 가장 적합한 계수를 산정하여 최종적으로 결정하였다. 각 경험계수 조건에 따른 선 공급원으로부터 1.5 m 높이에서의 중심선 속도의 오차 크기를 Fig. 14와 같이 나타내었다. 중심선 속도에 가장 민감하게 작용하는 경험계수는 TLEN, We_crit 및 Ca_crit로 나타났다.

Fig. 14.

Error of the centerline velocity at z = 1.50 m under each simulation conditions.

결과적으로 본 연구에서 수치 모델 해석을 위해 최종적으로 선정된 계수를 정리하면 항력계수(C_D)의 경우 기포입자를 구체로 가정하여 0.47의 값을, Richardson-Zaki 상관방정식의 상관 계수(C_rz)는 0.05, We_crit 및 Ca_crit의 경우 수중에서 방출되는 기포량에 따라 조정이 필요하다고 판단된다. 그리고 TLEN의 경우 computed 조건 또는 특정값(수심의 7% 이상)에서 흔들림이 발생하지 않았으며 본 연구에서는 유속 프로파일을 가장 잘 재현하는 computed 조건을 적용하였다.

상술한 수치해석 모델의 경험계수 민감도 분석 결과를 바탕으로 최적의 경험계수 값들을 적용하여 Kobus[1968]의 상승 기포로 유도된 중심선 속도의 계산 결과를 Fig. 15에 나타내었다. 경험계수들의 최적값은 Table 3과 같으며 수치해석 결과, 본 연구에서 적용된 공기 연행 및 드리프트 흐름 모델이 수중에서 방출되는 압축기포에 의한 흐름장 현상을 대체적으로 양호하게 재현할 수 있다고 판단된다.

Fig. 15.

Centerline velocity (wm) versus distance above air source (z) for different airflow rate (Qa).

Table 3.

Numerical model variables for sensitivity analysis

3.5 기포장막 주변의 유동장 해석

구축된 수치해석 모델을 이용하여 기포 방출량 조건에 따른 기포장막 주변의 유동장 재현 및 기포장막으로 유도되는 수평 유속 흐름의 크기를 계산하였다. 실험안은 Table 4와 같이 총 7가지의 조건으로 구분하였다. 수치해석 실험중 하나의 결과로서 Fig. 16은 공기방출량 372 LPM 조건에서의 기포장막 주변 유동장의 유속 크기 및 벡터를 나타내었다. 그림에서 상승하는 기포로 유도되는 기포장막 중심속도 및 표면에서 수평으로 편향하는 수평 유속 흐름의 재현이 잘 나타났으며 중심선을 기준으로 대칭으로 흐름장이 형성되는 것을 확인할 수 있다.

Table 4.

Numerical simulation cases

Fig. 16.

Simulation using the builted numerical model at Qa = 372 LPM.

또한 Fig. 17은 기포 방출량 180 LPM 조건에서 공급원에서의 거리 0.8, 1.1, 1.4 m 떨어진 지점에서의 밀도 및 수직 유속 프로파일을 도식화하였다. Table 4의 모든 실험안에서 기포 코어 너비와 기포장막 전체 폭의 비 λ는 0.58~0.63의 범위로 나타났으며 이는 Fannelrop and Sjeren[1980]과 Milgram[1983]가 제시한 범위와 일치하였다. 그리고 최대 수평 유속은 공급원 깊이 대비 0.15~0.23H(H는 수심) 지점에서 발생하였다. 또한 최대 수평 유속의 10%가 발생되는 거리를 외부 흐름 영역의 크기로 정의하면 공급원의 수면으로부터 깊이 대비 약 2.0~2.3H가 영향 범위로 판단된다.

Fig. 17.

Vertical velocity and density profile at Qa = 180 LPM.

본 연구에서 수행된 결과를 바탕으로 기포 방출량에 따른 최대 수평 유속은 선행 연구자들의 실험적 경험식과의 비교를 통해 Fig. 18에 나타내었다. 수치해석 시뮬레이션의 결과는 네덜란드 Delft 연구소에서 수행된 Abraham and Burgh[1962]의 실험 결과와 일치하는 것으로 판단된다.

Fig. 18.

Maximum horizontal velocity (um,max) for different airflow rate (Qa).

그리고 Fig. 19는 기포 방출량 변화에 따른 최대 수직 및 수평 유속의 관계를 나타내었는데 수치 모델을 이용한 결과를 이용하여 식(21)과 식(22)의 관계식을 산출하였다. 최대 수평 유속의 크기(u_m,max)가 최대 수직 유속 크기(w_m,max)보다 약 16.3%(=|w_m,max−u_m,max|/w_m,max) 작게 나타났는데 이는 연직 방향 흐름이 표면에서 수평 방향으로 전향될 때 에너지 손실에 따른 영향으로 판단된다. 이와 같은 결과는 전술한 Table 1의 Abraham and Burgh[1962]의 연구결과와 유사하다. 상술한 최대 수직 유속과 기포 방출량과의 상관성을 제시한 관계식은 향후 기포장막을 현장 적용함에 있어서 공학적인 기초 설계인자로 활용가능하리라 판단된다.

$\[ u_{m,\mathrm{m}\mathrm{a}\mathrm{x}}=1.164{\left(g{Q}_a\right)}^{1/3} \]$

(19)

$\[ w_{m,\mathrm{m}\mathrm{a}\mathrm{x}}=1.391{\left(g{Q}_a\right)}^{1/3} \]$

(20)

Fig. 19.

Maximum horizontal (um,max) and vertical (wm,max) velocities for different airflow rate (Qa).

4. 요약 및 결론

본 연구에서는 전산유체역학 상용프로그램인 FLOW-3D를 이용하여 수중 압축공기의 방출로 유도되는 기포장막의 수리학적 현상을 재현하기 위한 수치해석을 수행하였다. 수치모델의 적용에서 기포-물 상호작용과 난류 해석을 수행하고 기포 장막 해석을 위한 경험계수에 대해 민감도 분석을 수행하였다. 또한 최종적으로 기포공급량에 따른 기포장막의 거동과 주변 유동장에 대한 수리학적 특성을 분석하였다. 본 연구에서 얻어진 결과를 요약하면 다음과 같다.

(1) 본 연구에서는 FLOW-3D를 이용하여 수중 압축공기의 방출로 유도되는 기포장막의 모델 경험계수에 대해 민감도 분석을 수행하고 수리학적 현상을 재현하기 위한 수치해석모델을 구축하였다.

(2) 기포 연행 및 Drift-flux 모델을 적용하기 위하여 고려해야 하는 경험계수로는 총 6개로서 각 격자에 기포의 연행율을 결정하는 연행계수, 구체 형상으로 가정하는 기포의 항력계수, 초기 기포 크기, 드리프트 흐름에서 상대 속도 크기에 기여하는 Richardson-zaki 상관계수, 임계 웨버수 및 임계 카필러리 수에 해당하였다.

(3) 최종적으로 선정된 계수로 항력계수의 경우 기포입자를 구체로 가정하여 0.47의 값을, Richardson-Zaki 상관방정식의 상관계수는 0.05값으로 산정되었으며 We_crit 및 Ca_crit의 경우 수중에서 방출되는 기포량에 따라 조정이 필요함을 알 수 있었다.

(4) 방출된 기포와 유체 사이에 발생되는 난류 흐름을 계산하기 위하여 RNG k−ε 모델을 적용하였으며 이때 영향을 미치는 계수는 TLEN 이었다. TLEN 값이 작을수록 점성 영향으로 인해 상승 기포의 흔들림(wandering) 현상이 나타났다.

(5) 상술한 모델 경험계수를 바탕으로 기포 공급량에 따른 기포장막의 재현과 상승 기포로 유도되는 주변 유동장을 계산한 결과, 기포 코어 너비와 기포장막 전체 폭의 비 λ는 0.58∼0.63의 범위로 나타났으며 이는 선행연구 결과와 일치하였다.

(6) 기포장막 주변에서의 최대 수평 유속 발생 지점은 공급원 깊이(H) 대비해서 약 0.15∼0.23H 지점이며 최대 수평 유속의 10%가 발생되는 거리는 약 2.0∼2.3H 지점에 해당하였다.

(7) 수치해석 모델을 이용하여 기포 방출량 변화에 따른 최대 수평 및 수직 유속의 관계식을 산출하였는데, 최대 수평 유속은 u_m,max=1.164(gQ_a)^1/3, 최대 수직 유속은 w_m,max=1.391(gQ_a)^1/3의 관계식으로 제안할 수 있다. 최대 수평 유속의 크기가 최대 수직 유속 크기보다 약 16.3% 작게 나타났는데 이는 연직 방향 흐름이 표면에서 수평 방향으로 전향될 때 에너지 손실에 따른 영향으로 판단된다.

한편, 기포장막의 수리학적 현상을 재현하기 위한 수치해석모델은 향후 항만 및 항구 지역에 내습하는 파랑에너지 저감, 해양에서의 기름이나 부유물질 등의 오염물질 확산 저감, 발전소 취수구의 해양생물 퇴치 기술 등 다양한 연구 분야에서 기초 설계 모델로 활용가능할 것으로 판단된다. 하지만 실제 공학적인 설계에 적용할 수 있는 수준에 도달하기 위해서는 다양한 수심 및 기포 방출량 조건에서의 기포장막 주변 유동 해석에 관한 추가 연구가 필요할 것으로 생각된다.

Acknowledgments

이 논문은 2021년도 정부(해양수산부) 재원으로 영남씨그랜트센터의 지원을 받아 수행된 연구입니다(과제번호 : 2017046).

References

• Abraham, G. and Burgh, P., 1962, Reduction of salt water intrusion through locks by pneumatic barriers, Technical Report. Rijkswaterstaat and Delft Hydraulics Laboratory.
• Andreopoulos, J., Praturi, A., and Rodi, W., 1986, Experiments on vertical plane buoyant jets in shallow water, J. of Fluid Mechanics, 168(1), 305-336. [https://doi.org/10.1017/S0022112086000393]
• Batchelor, G.K., 1951, Note on a class of solutions of the Navier-Stokes equations representing steady rotationally-symmetric flow, Quart. J. of Mech. and Appl. Math., 4(1), 29-41. [https://doi.org/10.1093/qjmam/4.1.29]
• Brethour, J.M. and Hirt, C.W., 2009, Drift model for two-component flows, FLOW-3D Tech. Note FSI-09-TN83Rev., Flow Science, Inc., USA.
• Brevik, I., 1976, Partial wave damping in pneumatic breakwaters, J. of Hydraul. Div. 102(9), 1167-1176. [https://doi.org/10.1061/JYCEAJ.0004604]
• Brevik, I. and Kristiansen, O., 2002, The flow in and around air-bubble plumes, International J. of Multiphase Flow 28, 617-634. [https://doi.org/10.1016/S0301-9322(01)00077-5]
• Bulson, P.S., 1961, Current production by an air curtain in deep water, Dock Harbour Auth., 42, 15-22.
• Bulson, P.S., 1963, Large scale bubble breakwater experiments, Dock Harbour Auth., 44, 191-197.
• Bulson, P.S., 1968, The theory and design of bubble breakwaters, in proc. of Coast. Eng. Conf., ASCE, London, UK, II, 995-1015. [https://doi.org/10.9753/icce.v11.64]
• Ditmars, J.D. and Cederwall, K., 1974, Analysis of air-bubble plumes, in proc. of Coast. Eng.. Conf., Copenhagen, Denmark, 2209-2226. [https://doi.org/10.1061/9780872621138.131]
• Fannelop, T.K. and Sjoen, K., 1980, Hydrodynamics of underwater blowouts, AIAA 18th Aerospace Sciences Meeting, California, USA. [https://doi.org/10.2514/6.1980-219]
• Fannelop, T.K., Hirschberg, S. and Kuffer, J., 1991, Surface current and recirculating cells generated by bubble curtains and jets. J. Fluid Mech. 229, 629-657. [https://doi.org/10.1017/S0022112091003208]
• Isfahani, A.H.G. and Brethour, J.M., 2009, On the implementation of two equation turbulence models in Flow-3D, FLOW-3D Tech. Note FSI-09-TN86, Flow Science, Inc., USA.
• Jang, S.C., Kim, M.S., Lee, S.H. and Yoon, H.S., 2020, An experimental study on the effect of decreasing the surface wave energy due to the bubble curtain, J. Korean Soc. Mar. Environ. Energy, 23(4), 1-11. [https://doi.org/10.7846/JKOSMEE.2020.23.4.258]
• Kobus, H.E., 1968, Analysis of the flow induced by air-bubble systems, in proc. of Coast. Eng. Conf., London, UK, II, 1016-1031. [https://doi.org/10.9753/icce.v11.65]
• Kurihara, M., 1965, On the Study of a Pneumatic breakwater in Japan. Coast. Eng. in Japan, 8(1), 71-83. [https://doi.org/10.1080/05785634.1965.11924660]
• Levich, V.G. 1962, Physicochemical hydrodynamics, Prentice Hall, Englewood Cliffs, NJ, 464.
• Milgram, J.H., 1983, Mean flow in round bubble plumes. J. Fluid Mech. 133, 345-376. [https://doi.org/10.1017/S0022112083001950]
• Riess, I.R. and Fannelop, T.K., 1995, On deflecting drifting icebergs, in proc. of the 14th Int. Conf. Offshore Mechanics and Arctic Eng., Copenhagen.
• Schladow, S.G., 1992, Bubble plume dynamics in a stratified medium and the implications for water quality amelioration in lakes, Wat. Resour. Res., 28(2), 313-321. [https://doi.org/10.1029/91WR02499]
• Taylor, G.I., 1955. The action of a surface current used as a breakwater, in proc. of R. Soc. London, UK, Ser. A 231, 466-478. [https://doi.org/10.1098/rspa.1955.0188]
• Yakhot, V., Orszag, S.A., Thangam, S., Gatski, T.B. and Speziale, C.G., 1992, Development of turbulence models for shear flows by a double expansion technique, Phys. cs of Fluids 4(7), 1510-1520. [https://doi.org/10.1063/1.858424]