1. 서 론

공학 분야에서 높은 레이놀즈수를 갖는 난류 흐름 문제를 해석하기 위해서 일반적으로 사용되는 방법은 통계학적 난류 모형(statistical turbulence model)을 결합한 RANS (Reynolds-averaged Navier-Stokes) 모델링이다. 경계층(boundary layer)의 영향에 지배되는 난류 흐름에 대해서 계산의 정확도와 계산비용을 고려할 때, 계산상 가장 중요한 것은 경계층을 처리하는 방법이다(Durbin and Pettersson Reif, 2011). 높은 레이놀즈수 조건에서 벽면(solid wall) 근처의 흐름은 큰 유속 경사를 보이므로 마찰경사를 정확하게 계산하는 것은 어려운 문제이다. RANS 모델링을 할 때 벽면에서의 전단응력(shear stress)을 결정하는 방법은 두 가지가 있다. 첫째는 비활(no-slip) 조건을 적용하여 벽면에서 속도경사를 직접 결정하여 벽면에서의 전단응력을 결정하는 방법이다. 둘째는 대수층(log-layer)에 적용할 수 있는 반경험(semi-emperical) 공식을 이용해서 벽면에서의 전단응력을 계산하는 벽함수(wall functions)를 적용하는 방법이다.

난류모형(turbulence models)은 벽 근처에서의 난류흐름을 처리하는 방법에 따라 낮은 레이놀즈수(LRN, low Reynolds number) 모형과 높은 레이놀즈수(HRN, high Reynolds number)모형으로 나눌 수 있다. LRN 모형은 벽면에 접근할 때 다양한 흐름량의 제한된 거동(limiting behaviors)을 양호하게 재현하는 난류모형을 말한다. 여기서 LRN은 거시적 규모에서 흐름을 나타내는 것이 아니라, 점성 영향이 지배적인 벽면에 가까운 영역을 나타낸다. 즉, 점성저층부터 완충층(buffer layer)까지를 포함하는 전체 경계층을 세밀한 계산격자상에서 계산하여 경계층 흐름의 정확한 제한적 거동(correct limiting behavior)을 모의할 수 있는 것이 LRN 난류모형이다. 대부분의 난류운동에너지(turbulent kinetic energy)의 비소선율(specific dissipation rate) ω를 이용하는 모형들은 LRN 모형이다. 표준 $k-\epsilon$ 모형을 포함한 대부분의 일반적인 $k-\epsilon$ 모형들은 HRN 모형으로서 경계층을 직접 해석하는 대신 벽함수를 이용하여 대수층부터 계산을 수행한다. 하지만, 정확한 제한적 거동을 재현하는 소위 감쇄함수(damping functions)로 보완된 $k-\epsilon$ 모형은 LRN 모형에 포함된다. LRN 모형들은 경계층을 높은 정확도로 계산할 수 있지만, 벽면 근처에서 높은 해상도의 격자 구성이 필요하다. 이것은 높은 정확도를 위해서는 높은 계산비용이 필요함을 의미한다(Bredberg, 2000).

높은 해상도와 LRN 난류모형을 이용한 직접 계산을 통해서 벽면에서의 전단응력을 결정하는 방법이 계산의 정확도와 유체 물리학적 측면에서 선호되는 방법이지만, 벽 근처의 흐름을 해석할 때 벽면에 수직인 방향의 격자간격을 매우 작게(보통 무차원 벽거리 y⁺$\simeq$ 1) 을 유지를 해야 한다. 결과적으로 많은 계산격자가 필요하며, 계산시간 간격을 줄여야 하고, 계산격자의 가로-세로 비율을 너무 크게 하면 수렴 문제를 유발할 수 있다. 벽함수 경계조건을 이용할 경우는 벽면 근처 계산격자의 크기를 현저히 크게(보통 30 < y⁺ < 300 정도까지) 증가시킬 수 있어 전체 계산격자의 수를 현저히 줄일 수 있으며, 계산시간 간격을 크게 할 수 있어, 반복계산시 빠른 시간내에 수렴된 해를 구할 수 있다(Eça et al., 2015). HRN 조건에서 흐름분리(flow separation)와 재부착(reattachment)문제가 대상 흐름을 지배하지 않는 한 이러한 벽함수 방법을 적용하는 것이 일반적이다(Pope, 2000). 하지만 벽함수를 적용할 경우, 벽함수 경계조건은 점성흐름을 해석하기 위해서 적용할 수 있는 수학적 모형의 일부이며, 모델링의 정확도를 낮추는 대신 수치해석(numerical solution)을 쉽게 하는 것이 사용 목적이다.

벽함수 방법을 적용할 때는 난류모형에 따라 다양한 벽함수들의 조합을 이용할 수 있다. HRN 모델링을 위한 벽함수 연구는 다양한 유체역학 연구에서 많이 이뤄지고 있으며, 특히 선형 점성저층(viscous sub-layer)과 대수층(log-layer) 관계식을 혼합하는(blended) 방법이 출현하면서 벽함수의 사용이 증가하고 있다(Menter et al., 2003; Craft et al., 2006; Popovac and Hanjalic, 2007; Hedlund, 2014; Bäckar and Davidson, 2017; Lee and Paik, 2017).

이 연구에서는 수공학 문제에 널리 적용되고 있는 다양한 난류모형(turbulence models)들의 종속변수인 난류운동에너지(turbulent kinetic energy, TKE) k, 난류 소산율(turbulent dissipation rate) $\epsilon$, 난류에너지의 비소선율(specific dissipation rate) ω 그리고 와점성계수 ${\nu }_t$에 대한 벽함수들을 조합하여 각 기법의 성능을 평가하고자 한다. 모든 계산은 오픈소소 소프트웨어인 OpenFOAM 툴박스(OpenFOAM, 2018)를 이용하여 수행하였다.

RANS 수치모의의 정확도를 향상하기 위해서 벽면 근처에서 격자 간격을 세밀하게 유지하는 LRN 격자를 이용할 수도 있지만, 적용성이 떨어지므로, 이 연구에서는 실제 적용이 용이한 벽면 근처에서 균일 간격의 격자를 대상으로 분석을 수행한다. 먼저, 점성저층부터 대수층 상부까지 무차원 거리 y⁺ 값에 따른 무차원 유속 u⁺ 의 관계를 비교 평가한다. 그다음 후방계단(backward-facing step) 흐름 문제에 대한 모의 결과를 Durst and Schmitt (1985)의 실험값과 비교하여 모의 결과의 격자 해상도에 대한 민감도, 수치해의 수렴특성, 모형별 벽함수의 성능 특성을 제시한다.

2. 벽 근처 흐름영역과 유속분포

벽 근처에서의 평균유속분포를 무차원 유속과 무차원 거리의 항으로 그려보면, Fig. 1에서와 같이 여러 층(layers)과 영역(regions)으로 나눌 수 있다. 각 영역과 층은 벽으로부터의 거리 y를 경계층의 두께 $\delta$로 무차원화한(non-dimensionalized) 거리 그리고 y를 마찰유속(shear velocity) $u_{\tau }$와 동점성계수(kinematic viscosity) $\nu$로 무차원화한 벽단위(wall unit) y⁺ = $y{u}_{\tau }/\nu$의 함수로 구분된다. 마찰유속 $u_{\tau }$는 $\sqrt{{\tau }_w/\rho }$로 정의되며, 점성길이규모(viscous length scale)은 ${\delta }_{\mathit{\nu }}=\nu /u_{\tau }$이며, 마찰레이놀즈수(friction Reynolds number)는 $R{e}_{\mathit{\tau }}\mathit{=}{u}_{\mathit{\tau }}\delta /\nu =\delta /{\delta }_{\nu }$로 정의된다. 유속 u를 $u_{\tau }$로 무차원화한 u⁺ 를 y⁺ 에 대해서 선형적으로 그리면 Fig. 1과 같이 우리가 잘 알고 있는 벽근처의 일반적인 유속분포의 형태인데, y⁺ 축을 대수 스케일로 그리면 아래 그림과 같은 독특한 유속분포 그림이 된다.

Fig. 1.

Near-wall profiles of mean velocity plotted on (a) a linear scale and (b) a semi-logarithmic scale.
[Replotted from Kim et al., 1987 and Oweis et al., 2016]

벽면에서 유체 유속은 0이며, 벽으로부터 가장 가까이 위치한 층인 점성저층(viscous sublayers, y⁺ ≤ 5)에서는 난류의 영향은 무시할 수 있으며, 유속은 벽으로부터의 거리에 따라 선형으로 변한다. 점성저층과 대수영역(y⁺ > 30) 사이에 있는 영역을 완충층(buffer layer)이라고 하며, 이 층은 점성이 지배적인 흐름 부분과 난류가 지배적인 흐름 부분 사이의 천이영역(transition region)이다(Pope, 2000). 그림에서와 같이 대수법칙영역은 Re > 3,000 조건에서 30${\delta }_{\nu }$ < y < 0.3$\delta$범위에 존재한다. 레이놀즈수가 20,000을 초과하면 30${\delta }_{\nu }$ < y < 0.3$\delta$범위에 중첩층(overlap)이 있다(Perlin et al., 2016).

일반적인 벽근처에서의 유속분포는 점성영역(y⁺ < 5) 관계식과 대수영역(30 < y⁺ < 300) 관계식으로 나눌 수 있다. 점성영역은 유체가 점성의 영향을 받기 때문에 유체의 전단응력은 벽 전단응력 ${\tau }_w$과 같다고 볼 수 있다. 대수영역은 유체가 난류의 영향을 받고, 속도는 거리 y에 따라 log 함수로 매우 천천히 변화한다. 점성영역과 대수영역에 대한 벽함수는 각각 다음 식과 같이 정의된다.

여기서, $\kappa$는 von Kármán 상수(≒0.41)이며 B는 상수로서 매끄러운 벽면에 대해서 일반적으로 5.0을 적용한다.

3. 난류모형과 벽함수

3.1 난류 모형(turbulence model)

이 연구에서 사용한 난류모형은 공학적으로 널리 사용되고 있는 2방정식 모형인 표준 $k-\epsilon$모형(Launder and Sharma, 1974)과 난류 운동에너지와 비소산율이 종속변수인 k－ω SST(shear stress transport) 모형(Menter, 1993), Navier-Stokes 방정식을 다시 표준화하기 위해 RNG(Re-normalization Group) 방법을 사용하여 개발된 RNG $k-\epsilon$모형(Yakhot et al., 1992) 그리고 상대적으로 보다 최근에 개발된 1방정식 모형인 S-A (Spalart and Allmaras, 1992) 모형이다.

표준 $k-\epsilon$ 모형은 HRN 모형으로 난류 조건에 대한 평균 유동 특성을 모의하기 위해 전산 유체 역학에 사용되는 가장 일반적인 모델이다. $k-\epsilon$ 모델은 유동의 난류 특성을 나타내는 두 개의 종속변수인 난류 운동 에너지 k와 에너지소산율 $\epsilon$의 수송방정식을 계산하여 난류를 해석한다. 이 모델은 상대적으로 작은 압력 경사(gradient)를 갖는 자유 전단층 유동에 유용하다. 마찬가지로, 벽 경계 및 내부 흐름의 경우, 모델은 평균 압력 경사가 작은 경우에만 양호한 결과를 제공한다(Launder and Sharma, 1974).

k－ω 모형은 $k-\epsilon$ 모형과 마찬가지로 가장 보편적으로 사용되는 일반적인 2방정식 난류모형으로 두 개의 변수인 k와 ω에 대한 수송방정식을 풀어서 난류점성을 모의한다. k－ω모형은 LRN 모형이지만 벽 근처에서 해상도가 낮은 격자를 이용할 경우는 벽 함수와 함께 사용할 수도 있다. 이 연구에서 사용할 난류 모델은 Wilcox (1988)의 k－ω에서 발전된 모델인 SST (shear stress transport) 모형이다(Menter, 1993). k－ω SST 모델은 $k-\epsilon$모형과 k－ω 모형을 결합한 모델이다. 혼합 함수인 F₁은 벽 가까이에서 Wilcox 모형을 그리고 벽에서 떨어진 자유 흐름 영역에서 $k-\epsilon$ 모형을 활성화한다.

RNG $k-\epsilon$ 모형은 Yakhot et al. (1992)이 작은 규모(small scales) 유동의 효과를 구현할 목적으로 Navier-Stokes 방정식을 다시 표준화하기 위해 RNG (Re-Normalization Group) 이론이라 불리는 통계적 기법을 사용하여 개발하였다. 표준 $k-\epsilon$ 모델에서 와 점성(eddy viscosity)은 단일 난류길이규모(turbulent length scale)로부터 결정되므로, 계산된 난류확산은 설정된 규모에서만 발생한다. 하지만, 모든 규모의 유동들이 실제 난류확산에 기여하며, RNG 접근법은 난류 생성항에 대한 변화를 통해서 서로 다른 규모의 유동을 고려하는 수정된 형태의 $\epsilon$ 방정식을 이용한다. 아울러 표준 $k-\epsilon$ 모형은 난류 프랜틀(Prandtl) 수에 대해서 사용자가 설정하는 상숫값을 이용하는 한편 RNG 이론은 난류 프랜틀(Prandtl) 수에 대해서 해석적 공식을 제공한다. 표준 $k-\epsilon$ 모형은 HRN 모형이지만 RNG 모형은 LRN의 영향을 고려하는 유효점성에 대해서 해석적으로 유도된 미분식을 이용한다.

S-A 모형은 원래 공기역학 흐름을 해석하기 위해서 운동학적 난류점성에 대한 수송방정식을 해석하는 1 방정식 모형으로 역압력경사(adverse pressure gradient) 조건을 받는 경계층을 양호하게 해석한다. 지배방정식은 점성저층식과 대수층 관계식을 y⁺에 근거해서 혼합한 형태를 보인다. 따라서 점성저층 모의를 위해서 추가적인 비선형 감쇄함수를 이용하지 않는다. 본래 이 모형은 실질적으로 LRN 모형이며, 경계층의 점성에 영향을 받는 영역을 적절하게 해결해야 한다. 점성저층 관계식에서 y⁺에 기초한 로그 공식으로 자동 혼합(blending)한다. 이러한 혼합은 완충층(5 < y⁺ < 30)에서 중간 정도의 값을 적절히 처리하도록 검증되었다. 벽면에 부착된(attached) 흐름과 완만한 흐름 분리(flow separation) 및 재순환(recirculation)을 포함한 흐름 해석에 경제적이면서 정확하지만 강한 흐름분리 및 자유전단흐름 해석에 단점이 있는 것으로 알려져 있다.

3.2 벽함수(wall functions)

난류모형들에서 종속 변수인 k, $\epsilon$, ω 그리고 난류점성계수 ${\nu }_t$에 대해서 적용 가능한 벽함수들의 성능을 평가하고자 한다. 적용한 벽함수들을 간략하게 요약하면 아래와 같다. 추가적인 설명을 위해서는 Liu (2016)을 참고할 수 있다. 이 절에서 사용하는 각 벽함수의 명칭은 약칭은 사용하였으며, 실제 OpenFOAM toolbox에서 사용하는 명칭은 괄호안에 명기하였다.

3.2.1 k-벽함수

난류운동에너지 k에 대한 경계값 설정은 난류 모델링 방법에 따라 다음의 두 가지 벽함수를 이용할 수 있다. 이하 벽함수 명칭 중 괄호 안은 OpenFOAM에서 해당 벽함수의 명칭을 나타낸다.

kq-벽함수(kqRWallFunction)

HRN 모델링에서 k에 대한 경계조건은 아래의 관계식에 근거한 벽함수를 이용하여 설정한다.

여기서 $C_{\mu }$=0.09이다.

kL-벽함수(kLowReWallFunction)

LRN 모델링에서도 적용이 가능한 벽함수이며, 기본적으로 계산에서 $y_L^{+}(=\log (\max \left[E\times C,1\right])/\kappa )$ 값을 기준으로 점성 및 대수 영역을 구분하여 다음과 같이 계산한다.

점성 영역 $y^{+}<y_L^{+}$

여기서, $k^{+}=k/u_{\tau }^2$이며 C_eps2 = 1.9, $C_f=1/{(y++C)}^2+2y^{+}/c^3-1/c^2$, C = 11.0이다.

대수 영역 $y^{+}>y_L^{+}$

여기서, C_k = -0.416, B_k = 8.366이다.

3.2.2 $\epsilon$-벽함수(EpsilonWallFunction)

k 벽함수와 다르게 $\epsilon$ 벽함수는 계산격자 면(face)에서의 값이 아닌 계산격자 중심에서의 값을 계산한다. 여기서 사용되는 매개변수 W는 하나의 셀에서 몇 개의 면(faces)이 경계조건 설정에 이용되는지를 나타낸다. 예를 들면 2차원 계산에서 사각형 계산격자를 이용한다면 셀의 흐름방향 상하류 2개 면과 셀의 상부 면을 포함한 3개 면에 대해서 벽함수를 적용하므로 W는 3이고, 3차원 계산의 경우 6면체 계산격자를 이용한다면 바닥 벽면을 제외한 나머지 5개 면에 대해서 벽함수를 적용하므로 W는 5가 된다. 계산격자의 각 면에서 $\epsilon$값은 y⁺ 값에 따라서 다음과 같이 계산된다.

점성 영역 $y^{+}<y_L^{+}$

대수 영역 $y^{+}>y_L^{+}$

여기서 난류운동에너지 k는 계산격자 중심에서의 값이고 y_i는 각 격자면의 중심에서 벽까지의 거리이다.

3.2.3 ω-벽함수(OmegaWallFunction)

ω-벽함수는 $\epsilon$-벽함수와 마찬가지로 셀 면에서의 값을 계산한 후 그것의 혼합으로 셀 중심의 값을 계산한다. 다음과 같은 점성저층과 대수층 방정식의 조합으로서 이루어진다.

점성 영역 $y^{+}<y_L^{+}$

대수 영역 $y^{+}>y_L^{+}$

여기서 ${\beta }_1$= 0.075이다. 대수 영역과 점성 영역의 중간 지역에서는 두 관계식을 다음과 같이 조합하여 ω값을 계산한다.

3.2.4 ${\nu }_t$-벽함수

난류점성계수 ${\nu }_t$의 경계값을 설정하기 위해서는 먼저 다음과 같이 표현되는 바닥전단응력을 고려해야 한다.

만약 DNS (direct numerical simulation)를 수행하여 해를 구했다면 바닥에서의 속도경사를 알 수 있으므로 Eq. (11)의 응력을 직접 계산할 수 있다. 그렇지 않으면 근사적으로 다음과 같이 구할 수도 있다.

여기서 u_c와 u_w는 각각 첫번째 계산격자의 중심 그리고 바닥면에서의 유속이다. Eq. (12)가 성립하기 위해서는 첫번째 계산격자의 중심이 난류의 영향을 무시할 수 있는 점성저층에 위치해야 한다. 하지만, 벽 근처에서 상대적으로 큰 격자간격을 적용하는 HRN RANS 모델링은 이러한 격자조건을 맞추기가 어려우므로 난류의 영향을 반영하도록 새로운 유효점성(effective viscosity) ${\nu }_{eff}$을 고려하여 바닥전단응력을 수정하여야 한다.

대수영역에서 무차원 난류운동에너지는 다음과 같이 표현된다.

이 식으로부터 마찰유속은 $u_{\tau }=C_{\mu }^{1/4}{k}^{1/2}$로 표현되며, 마찰유속은 다시 $u_{\tau }=u/(u/u_{\tau })=u/u^{+}$의 관계로부터 바닥전단응력 ${\tau }_w$(= $\rho u_{\tau }^2$)는 다음과 같이 표현할 수 있다.

Eq. (12)와 Eq. (14)의 관계로부터 ${\nu }_{eff}$를 다음과 같이 정의할 수 있다.

결과적으로 난류점성 ${\nu }_t$를 얻을 수 있다.

OpenFOAM에서는 기본적으로 Eq. (16)을 이용하여 의 경계조건을 설정하는데, 여기서 y⁺를 계산하는 방법 그리고 $u_{\tau }$를 정의하는 방법에 따라 유도된 다양함 벽함수들을 적용할 수 있다. 이 연구에서는 공학적으로 적용이 적합한 다음 4가지의 난류점성 벽함수들의 성능을 평가하고자 한다.

nutk-벽함수(nutkWallFunction)

nutk-벽함수는 대수영역($y^{+}>y_L^{+}$)에서 $u_{\tau }=C_{\mu }^{1/4}{k}^{1/2}$의 관계식에 의해서 유도된 Eq. (16)을 이용하여 값을 설정한다. 첫번째 계산격자의 중심이 점성영역($y^{+}<y_L^{+}$)에 위치하면 ${\nu }_t$ = 0으로 설정한다. 이 벽 함수는 벽 근처 첫번째 계산격자의 중심이 점성저층 또는 대수영역에서 있을 때만 효과적으로 계산되기 때문에 계산 오차를 줄이기 위해서는 첫번째 계산격자의 중심이 완충층에 위치하지 않도록 격자를 구성할 필요가 있다.

nutU-벽함수(nutUWallFuction)

nutU-벽함수는 무차원 거리 y⁺를 계산하는 방법을 제외하면, nutk-벽함수와 같은 방법으로 ${\nu }_t$ 경계조건을 설정한다. 대수영역에서 벽으로부터의 거리와 유속분포관계는 다음 같이 주어진다.

또는 다른 형식으로 다음과 같이 주어진다.

여기서 일반적으로 사용되는 경험상수 E와 B의 값은 각각 7.76과 5.0이다.

이 식은 다음과 같이 변형시킬 수 있다.

결국, 다음과 같이 정리된다.

Newton-Rapshon 반복법으로 Eq. (20)을 계산하여 y⁺ 을 구한 후 nutk-벽함수와 같게 ${\nu }_t$값을 계산한다. nutU-벽함수도 대수법칙 관계식을 적용하므로 nutk-벽함수를 이용할 때와 같이 첫 번째 계산격자의 중심이 대수영역에 위치하도록 격자를 구성할 필요가 있다.

nutUB-벽함수(nutUBlendedWallFunction)

이 벽함수는 점성저층에서의 마찰유속관계식과 대수영역에서의 마찰유속관계식을 혼합하여 구한 합성 마찰유속을 이용하여 Eq. (16)을 계산함으로써 난류점성의 경계값을 설정하는 방법이다(Menter et al., 2003).

여기서, 지수 n은 일반적으로 4를 적용하고 $u_{\tau ,vis}$와 $u_{\tau ,\log }$는 각각 점성저층과 대수영역에서의 마찰유속으로 다음과 같이 정의된다.

이 처리방법은 점성저층과 대수영역에서 각각 유효한 마찰유속을 혼합하여 적용하므로 LRN 모델링을 위한 점성관계식으로부터 HRN 모델링에서의 벽함수 처리방법으로 격자의 밀도에 따라 자동으로 변환된다. 따라서 첫번째 계산격자의 위치에 민감하지 않아, 완충층에서 발생하는 기존 벽함수의 성능 저하를 방지할 수 있다.

Spalding 관계식

Spalding (1961)이 제안된 y⁺ 와 u⁺ 의 관계식을 기반으로 하며, y⁺ 는 Newton-Raphson방법을 이용하여 계산한다. 결과적인 관계곡선식은 점성저층에서 u⁺ = y⁺ 그리고 대수영역에서는 $u^{+}=E{y}^{+}/\kappa$ 곡선에 맞춰지도록 적합된(fitted) 관계식을 적용한다.

여기서 $\lambda =\kappa u^{+}$이며, 연속적으로 변하는 y⁺ 와 u⁺ 의 관계를 사용하므로 첫 번째 격자의 중심이 어느 영역에 속하든지 상관없이 적용할 수 있다. 이 벽함수는 LES (large eddy simulation) 또는 LRN RANS 모델링을 수행할 때 지역적으로 y⁺ 값이 1보다 크면 발생하는 경계조건 오차를 줄이기 위해서 적용하는 것이 일반적이다.

4. 벽함수 성능 평가

난류모형별로 적용한 벽함수의 성능을 비교하기 위해서 다음의 두 단계 절차를 수행하여 무차원 유속과 무차원 거리의 관계를 도시하였다. 먼저 벽에 수직인 방향으로 40개의 계산격자를 이용하여 격자를 생성하고, 주어진 y (첫번째 셀 중심에서 벽면까지의 거리)에 대해서 유체의 점성계수를 변화시키고 계산을 수행하여 다양한 y⁺ 값에 대해서 유속경계값을 구한다. 모든 계산격자에서 유속의 초기조건은 단면평균 유속으로 설정하고 계산을 수행하였다. 따라서 벽 근처에서의 유속은 평균유속에서 점차로 줄어들면서 수렴하게 된다. 이렇게 벽 경계 계산격자에서 계산된 유속을 마찰유속으로 무차원화한 유속인 u⁺ 값들을 y⁺ 값에 대해서 Fig. 2와 같이 도시하였다. 그림에서 각 모델링을 통해 구한 값들을 Spalding 식으로 구한 곡선과 비교하여 평가하였다. 이 절에서는 벽함수를 적용한 다양한 난류모형을 이용한 모델링시 계산 결과의 수렴특성을 분석한 후, 난류모형별 적용한 벽함수의 계산 결과를 평가하고자 한다.

Fig. 2.

Velocity profiles computed by RNG k-ε model with nutk-WallFunction at difference modeling times and convergence criteria

4.1 해의 진동 특성

수치모델링의 전체 계산시간은 2000 s와 3000 s를 적용하였으며, 반복계산을 수행할 때 각 계산단계에서 반복해의 오차허용범위(tolerance)를 난류량과 유속의 잔차(residual)에 대해서 1×10^－6, 1×10^－9, 1×10^－12값을 적용하였다.

분석결과 가장 주목해야 하는 것은 RNG $k-\epsilon$모형을 nutk-벽함수와 함께 적용할 때, 첫번째 격자가 완충층에 놓이면 오차허용범위를 크게(예를 들면 1×10^－6) 설정할 경우 해가 진동을 하면서 비합리적인 경계값을 도출한다는 것이다. Fig. 2에서 보인 바와 같이, 해가 보다 더 수렴되도록 계산시간을 2000 s에서 3000 s로 증가시키더라도 각 계산단계에서 수렴조건을 충분히 줄여주지 않으면, 점성저층과 대수층 사이의 완층층에서는 해가 수렴하지 않고 진동한다. 이 경우 각 계산단계에서의 오차허용 범위를 1×10^－9 이하로 낮추거나 반복계산법 적용시 최소 반복계산 횟수를 3회 이상으로 설정하여야 안정적으로 수렴된 해를 구할 수 있는 것으로 나타났다. 한편, 표준 $k-\epsilon$모형과 k－ω SST 모형은 nutk-벽함수와 함께 적용하더라도 완층층에서 경계값이 진동하는 현상은 없는 것으로 나타났다. 그리고 RNG $k-\epsilon$ 모형도 nutU-벽함수와 함께 적용하는 경우 해의 진동현상은 없는 것으로 나타났다.

4.2 해의 수렴특성

완충층에서 해의 진동이 발생하지 않는 조건, 즉 잔차의 허용값을 1×10^－9로 설정하고 계산시간에 따른 해의 수렴특성을 비교하였다. 총 1500 s 까지 계산된 유속을 시간에 따라 도시하면 Fig. 3과 같다. HRN 난류모형의 경우 상대적으로 큰 오차가 발생하는 완충층 영역에서의 해의 거동을 보기 위해서 y⁺ 가 완전한 완충층에 속하는 10과 20, 대수영역의 경계에 있는 30 그리고 완전한 대수영역인 70에 대해서 계산을 수행하였다. 오차허용조건을 충분히 작게 설정한다면, $k-\epsilon$ 모형은 계산시간별 유속의 변화는 각 난류모델 별로 nutk-벽함수나 nutU-벽함수 둘 중 어는 것을 사용하여도 수렴특성은 큰 차이가 없으므로 nutU-벽함수를 적용한 결과만을 Fig. 3에 도시하였다. 대표적인 HRN 모형인 $k-\epsilon$ 모형은 첫번째 계산격자의 중심이 완충층에 위치하면 초기에 해가 시간에 따라 크게 진동한 후 안정화 되는 것으로 나타났다. SST 모형은 적용할 경우 계산 초기에 바닥 부근에서의 유속이 크게 한번 진동한 후 안정화되는 것으로 나타났다. S-A 모형은 이러한 진동 현상이 거의 발생하지 않는 것으로 나타났다.

Fig. 3.

Time series of velocity components computed by turbulence models at different four wall distances off the wall

수렴된 값의 특징을 보면 완충층에서는 $k-\epsilon$과 RNG $k-\epsilon$ 모형은 경계 유속이 상대적으로 크게 계산되었으며, SST 모형이 가장 작은 값으로 계산되었다. 한편, 첫번째 계산격자가 대수영역에 위치할 경우 RNG $k-\epsilon$과 SST는 매우 유사한 유속이 계산되었으며, $k-\epsilon$과 S-A보다 상대적으로 유속값은 작은 것으로 나타났다.

4.3 벽함수 해석 결과

수렴된 해를 이용하여 난류모형별 u⁺－y⁺ 관계를 도시하면 Fig. 4와 같다. 먼저 Fig. 4(a) and (b)를 비교해보면 $k-\epsilon$모형은 k에 대해서 kq-벽함수 또는 kL-벽함수 무엇을 적용해도 그리고 ${\nu }_t$에 대해서 nutk-벽함수나 nutU-벽함수 무엇을 적용해도 결과에는 차이가 없는 것을 알 수 있다. 반면, SST 모형은 nutk-벽함수를 적용하면 nutU-벽함수에 비해 완충층(y⁺ = 10~20) 영역에서 유속을 과대 산정하는 것으로 나타났다. RNG $k-\epsilon$은 Fig. 4(c)에서 보인 바와 같이, 충분한 수렴이 이뤄진다면 벽 근처 유속 값은 벽함수의 종류에 민감하지 않은 것으로 나타났다. Fig. 4(b) and (d)를 비교하면, SST 모형은 nutU-벽함수를 사용하는 것보다 Menter et al. (2003)이 제시한 혼합벽함수인 nutUB-벽함수를 이용하는 것이 완충층에서의 유속의 과대 산정을 억제하고, Spalding곡선에 근접하면서 점성저층과 대수층에서 상대적으로 부드럽게 천이되는 것으로 나타났다. 아울러 S-A모델도 SST 모형과 같이 nutUB-벽함수를 이용하면 보다 합리적인 유속값을 계산하는 것으로 나타났다.

Fig. 4.

Near wall velocity profiles computed by selected turbulence models with wall functions

5. 적 용

공학적으로 모델 검증에 널리 이용되고 있는 후방계단흐름(backward-facing step flow)에 난류모형들을 적용하여 성능을 평가하였다. 대상 흐름은 계단 높이 h와 평균유속에 근거한 레이놀즈수 Re = 155,000이며 계산된 재순환(recirculation) 흐름의 재부착(reattachment) 거리 x_r과 흐름 방향 유속 분포를 Durst and Schmitt (1985)의 실험값과 비교하여 모델링 결과의 특성을 분석하였다. 격자 민감도 분석을 위해 3가지 해상도의 격자, 즉 성긴격자(coarse mesh, 2416+5632), 중간격자(medium mesh, 4832+11264) 그리고 조밀격자(fine mesh, 7248+16896)를 사용하였다. 이들 중 성긴격자를 도시하면 Fig. 5와 같다.

Fig. 5.

Computational mesh for backward-facing step flow

유입부에서 완전히 발달한 흐름을 재현하기 위하여 매핑(mapping) 기법을 사용하였다. 매핑기법은 유입 단면에서 하류로 $2\pi h$이상 떨어진 단면에서 계산된 수리량과 난류량들(u, k, $\epsilon$, ω, $\tilde{\nu }$, 등)을 매 계산단계마다 반복적으로 복사하여 유입 경계조건으로 사용하여 완전 발달한 흐름을 재현하는 방법이다(Kaikstis et al., 1991). 단면이 일정한 직선 수로에서 적용할 수 있는 기법으로서 이 연구에서는 입구부로부터 8h (여기서 h는 계단 높이)만큼 떨어진 하류 단면을 이용하여 매핑기법을 적용하였다. 이 절에서 제시된 결과는 k에 대해서는 모두 kL-벽함수를 적용한 결과이다. 난류점성계수 ${\nu }_t$에 대한 벽함수는 $k-\epsilon$, RNG $k-\epsilon$, k－ω SST 모형은 nutk-벽함수를 그리고 SST-nutUB와 S-A 모형에 대해서는 nutUB-벽함수를 적용한 결과이다.

해상도가 다른 3개의 격자에서 계산된 수치모의 결과를 토대로 벽면의 첫번째 계산격자에서 계산된 무차원 거리의 최소, 최대 그리고 평균값을 정리하면 Table 1과 같다. 당연한 결과이지만, 격자의 해상도를 높일수록 y⁺ 값은 감소함을 확인할 수 있다. 평균 y⁺ 값은 성긴격자에서 100~200범위의 값을 가지며, 조밀격자에서는 33~57 사이의 값으로, 모두 대수영역에 포함된다고 볼 수 있다. 재순환 흐름과 관련된 재부착 지점에서는 벽면에 평형한 유속이 매우 작은 값을 가지므로, 이 부근에서 계산된 최소 y⁺ 값은 조밀격자에서 모두 점성저층(y⁺ < 5)에 포함되는 것으로 나타났다. 적용한 난류모형별 특징을 보면, $k-\epsilon$과 RNG $k-\epsilon$ 모형은 SST-nutUB와 S-A 모형에 비해 y⁺ 가 큰 것으로 나타났다. 이것은 $k-\epsilon$ 계열의 모형이 벽 경계에서 상대적으로 큰 마찰유속 값을 계산하기 때문이다. SST 모형도 Menter et al. (2003)이 제시한 nutUB-벽함수가 아닌, nutk-벽함수를 적용할 경우 $k-\epsilon$ 모형과 유사한 y⁺ 값이 계산되는 것으로 나타났다.

Table 1. Dimensionless distances at the first cell center off the wall computed by applied turbulence models

대표 계산 결과로 두 개의 난류모형에 의해서 계산된 유속벡터장과 흐름방향 유속성분을 구간 -3 ≤ x ≤ 18에 대해서 도시하면 Fig. 6과 같다. 두 수치해 모두 계단 상류에서는 완전 발달한 전형적인 유속분포를 보이며 계단 모서리에서 발생하는 흐름분리에 의해서 생성되는 재순환 영역을 벡터장에서 볼 수 있다. 계단모서리에서 분리된 유체가 바닥면에 부착되는 재부착점(reattachment point)과 계단사이의 거리인 재부착거리 x_r은 실험에서 8.5로 관측되었다(Durst and Schmitt, 1985).

Fig. 6.

Velocity vectors and streamwise velocity component computed by [upper] the standard $k-\epsilon$ with nutk-wall functions and SST model with nutUB-wall functions on the coarse mesh

적용한 난류모형에 의해서 계산된 재부착거리를 정리하면 Table 2와 같다. 중간격자에서 계산된 값은 성긴격자와 조밀격자의 대략 중간값을 갖기 때문에 표에 포함하지 않았다. 표에서 주어진 바와 같이, 전반적으로 성긴격자보다는 조밀격자에서 계산된 결과가 Durst and Schmitt (1985)의 실험값과 잘 일치하는 것으로 나타났다. 표준 $k-\epsilon$ 모형을 제외한 나머지 모형들은 전반적으로 양호하게 재부착거리를 잘 예측하는 것으로 나타났다. 이러한 결과는 이미 잘 알려진 것으로서, 표준 $k-\epsilon$ 모형은 HRN 모형으로 흐름분리를 동반하거나 압력경사가 현저한 흐름은 적절히 해석할 수 없기 때문이다(Launder and Sharma, 1974). RNG $k-\epsilon$ 모형은 재부착거리를 양호하게 예측하며, 계산결과는 격자해상도에 민감하지 않은 것으로 나타났다. SST-nutUB 모의 결과는 조밀한 격자에서 실험값과 가장 가까운 값을 보여주는 한편, 해석결과는 다소 격자해상도에 민감한 것으로 나타났다.

Table 2. Reattachment lengths computed by applied turbulence models

선택된 4개의 단면에서 계산된 유속분포를 실험값과 비교하면 Fig. 7과 같다. Fig. 7(a)는 계단으로부터 2.5h 상류에 있는 단면으로 매핑기법에 의해서 완전 발달한 흐름이 재현된 지점으로 계산결과는 실험값과 잘 일치하고 있다. 이 결과는 매핑기법이 성공적으로 작용했음을 보여준다. Fig. 7(b)는 계단으로부터 0.24x_r 하류에 있는 지점으로 재순환 흐름 영역의 상류 부분이며 수로의 아랫부분에서는 역방향 흐름이 발생하는 것이 특징이다. 전반적으로 S-A 모형을 제외한 대부분 모형이 유사한 정확도로 실험값을 재현하고 있다. S-A 모형은 바닥 부근에서 역방향 흐름을 과대 산정하는 것으로 나타났다. Fig. 7(c)는 재부착 지점 바로 직상류 지점인 0.94x_r에서의 유속분포이다. 마찬가지로 S-A 모형을 제외한 다른 모형들은 유사한 계산결과를 보이며, 수로 상부 부분에서 실험값보다 다소 큰 유속을 발생시키고 있다. 이러한 최대유속의 과대산정은 Friedrich and Arnal (1993)의 LES (large-eddy simulation) 수치해와 유사하다. 재부착거리의 2.35배 하류 지점에서의 유속분포를 보여주는 Fig. 7(d)는 모든 수치해석 결과가 수로의 상부에서는 유속을 과대산정하고 수로의 하부에서는 유속을 과소산정하는 것으로 나타났다. 재순환 흐름의 강한 비정상성(unsteadiness) 때문에 재부착 거리 하류에서는 흐름이 빠르게 회복되어 유속분포가 수로 중앙을 중심으로 대칭을 이루게 되는데, 본 연구에서 수행한 RANS 수치모의에서는 이러한 흐름 회복을 합리적으로 재현하지는 못하는 것으로 나타났다. 추후 LES 또는 DES (detached-eddy simulation) 수치모의를 수행하여 이러한 흐름 회복의 특성을 연구하고자 한다.

Fig. 7.

Comparison of computed velocity profiles with measurements at selected 4 locations. In legend, kE = $k-\epsilon$; kO = $k-\omega$ SST; nk = nutk-wall function; UB = nutUBlended; SA = S-A; c = coarse mesh; f = fine mesh

6. 결 론

난류 흐름 계산을 위한 RANS 수치모의에서 널리 이용되고 있는 난류모형들은 벽함수와 함께 적용된다. 가용한 벽함수들을 이용하여 난류모형과 벽함수들의 성능을 비교·분석하였으며, 그 결과를 요약하면 다음과 같다.

RNG $k-\epsilon$ 모형을 nutk-벽함수와 함께 적용할 경우, 벽면으로부터 첫번째 계산격자의 중심이 완충층(6 < y⁺ < 11)에 위치하게 되면 수치적으로 안정된 경계값을 구하기가 어렵고 경계 유속이 과소산정될 수 있다. 이러한 오류를 제거하기 위해서는 반복계산시 허용오차(잔차의 허용값)을 1×10^－9 정도로 작은 값으로 설정하거나, 반복법을 적용할 때 최소 반복계산 수를 3회 이상으로 설정할 필요가 있다. 한편, RNG $k-\epsilon$ 모형도 nutU-벽함수와 함께 적용하면 이러한 오류는 없는 것으로 나타났다.

벽함수 해의 수렴특성은 nutk-벽함수와 nutU-벽함수 모두 유사한 것으로 나타났다. 표준 $k-\epsilon$ 모형은 첫번째 계산격자가 완충층(10 < y⁺ < 20)에 놓이게 되면 계산 초기 단계에 해가 크게 진동하는 것으로 나타났다. 이러한 결과는 일반적으로 첫 번째 계산격자의 중심이 대수층(또는 점성저층)에 위치하는 경우에 유효하도록 벽함수가 개발되었기 때문으로 판단된다. SST 모형을 적용한 경우 바닥 부근에서의 유속이 계산 초기에 한 번 크게 진동한 후 안정화 되는 것으로 나타났다.

표준 $k-\epsilon$이나 RNG $k-\epsilon$모형은 k대해서 kq-벽함수나 kL-벽함수 모두 그리고 ${\nu }_t$에 대해서는 nutk-벽함수나 nutU-벽함수 모두에 민감하지 않은 것으로 나타났다. SST 모형은 nutk-벽함수를 이용할 경우 첫번째 계산격자 중심이 완충층에 놓이면 다른 벽함수들에 비해 유속을 과대산정하는 것으로 나타났다. 한편, nutU-벽함수를 이용할 경우는 $k-\epsilon$이나 RNG $k-\epsilon$ 모형과 유사한 결과를 도출하는 것으로 나타났다. SST 모형을 적용할 때는 k에 대해서는 kL-벽함수를 그리고 ${\nu }_t$에 대해서는 nutUB-벽함수를 적용하는 것이 가장 우수한 결과를 도출하는 것으로 나타났다.

후방계단흐름에 적용한 결과, 표준 $k-\epsilon$ 모형은 약 13% 정도 재부착거리를 과소산정하는 것으로 나타났으며, 다른 모형들은 격자 해상도를 적절히 유지해준다면 약 3% 오차범위 안에서 재부착거리를 양호하게 예측하는 것으로 나타났다. RNG $k-\epsilon$ 모형은 재부착거리 측면에서 보면 계산결과는 격자해상도에 가장 민감하지 않은 것으로 나타났다.

후방계단흐름 수치모의 결과를 실험값과 비교해보면, 재순환 흐름 영역 상부에서의 큰 유속분포는 모두 양호하게 예측하는 반면에 바닥 부근에서 S-A 모형은 역방향 유속을 과대산정하고, RNG 모형 결과가 가장 실험값과 유사한 것으로 나타났다. 재부착지점에서의 유속분포는 $k-\epsilon$과 RNG $k-\epsilon$ 모형이 최대유속을 과대산정하는 한편, SST 모형이 기존 LES 해석결과와 가장 유사하면서 실험값을 양호하게 예측하는 것으로 나타났다. 재순환 영역 하류에서는 유속분포가 계단의 영향에서 벗어나 회복되면서 수로 중심을 기준으로 대칭을 이루는 흐름 회복이 특징인데, 모든 난류모형이 이러한 흐름 회복을 과소평가하는 것으로 나타났으며, 이것은 RANS 수치모의의 한계로 판단된다.