발간 현황

Chorin의 가상 압축성 기법[15]은 타원(eliptic) 형태를 가지는 연속방정식에 가상 시간 항을 추가함으로써 그 형태를 쌍곡(hyperbolic)형 방정식으로 변화시키므로, 비압축성과 압축성 간의 코드 전환이 용이하다는 장점을 가진다. Chorin의 가상 압축성 기법으로 표현된 비압축성 지배 방정식은 직교좌표계의 질량, 운동량 보존 방정식으로, 다음과 같이 표현된다.

여기서 β는 압축성 계수이고, τ는 가상 시간(Pseudo time)으로써 모두 수렴속도에 영향을 미치는 인자이다. 압력 p와 속도항들 u_i 모두가 가상 시간에 따라 변화가 없으면 Eq. 1과 Eq. 2는 원래의 보존 방정식 형태로 돌아오는 것을 알 수 있다.

액체-기체 경계면 추적에는 CLSVOF 방법을 사용하며, 경계면에서의 상변화 및 확산 과정을 무시하면 VOF, LS 이송 방정식은 각각 다음과 같다.

여기서 F는 VOF 값이고, ϕ는 LS 값이다. 이송 방정식으로부터 구한 체적 분율을 이용해 셀 내의 밀도 값, ρ = (1 - F) ρ_g + F_ρl 을 구한다.

Eq. 3의 차분화 (Discretization)는 다음과 같다.

VOF 값은 경계면 근처에서 불연속하기 때문에 Eq. 5, 6, 7의 Flux 값을 계산해 주기 위해선 Reconstruction 과정을 거쳐야 한다. 본 연구에서는 Reconstruction 방법 중 하나인 PLIC 방법을 적용하였다[16,17].

VOF 이송 방정식을 계산하는 과정에서 생기는 과도(Overshoot)과 과소(Undershoot) 오차는 다음과 같은 제한자로 조정한다.

Eq. 4는 LS 값이 연속적이기 때문에 다음과 같이 쉽게 차분화 할 수 있다.

Eq. 9, 10, 11은 2차 정확도 중심 차분법으로 계산하였다.

LS 이송 방정식 계산 뒤에, LS 값은 경계면에 대한 거리와 다르기 때문에 LS 값을 재초기화(Reinitialization) 과정이 필요하다. CLSVOF에서는 이런 재초기화 과정이 LS와 VOF를 연결시키는 과정이라 할 수 있다. LS 재초기화는 다음과 같은 거리 함수로 표현된 방정식을 반복계산을 통해 구할 수 있다.

여기서 ϕ₀은 초기 LS 값이고, τ는 재초기화에 대한 가상시간이고, ∇ϕ는 2^nd 중심 차분법으로 계산한다. sign(ϕ₀)는 다음과 같다.

Eq. 13을 계산하기 위해서는 경계면에서의 LS 값을 고정시켜야 하는데, 여기서 경계면의 LS 값은 VOF 값을 이용해 계산한다 [10]. 또한 계산 효율성을 위해 LS와 VOF 이송 및 재초기화 과정은 경계면이 있는 셀에서의 거리가 ±이내에 있는 셀에서만 계산한다.

Eq. 1, 2의 공간 차분화는 상류 차분법으로 계산하고, 상류 차분법의 Reconstruction 과정은 WENO[18], 그리고 Riemann 해법은 LDFSS[19]를 적용하였다. 시간 차분화는 수렴 가속화를 위해 LU-SGS (Lower-upper Symmetric Gauss-Seidel) 기법을 적용하였다.

개발된 CLSVOF 프로그램의 계산 검증을 위해 2차원 Zalesak’s disk[20]와 액체의 3차원 변형(Deformation) 계산을 실시했다. 계산 격자는 정육면체의 정렬격자를 사용했고, 격자 의존도를 파악하기 위해 격자의 길이(h)를 h = 1/64, 1/128, 1/160의 3가지에 대해서 계산을 실시하였다. 계산 효율성을 증대시키기 위해 MPI(Message Passing Interface) Parallel 계산 기법을 적용했고, 각 블록당 h = 1/32의 격자를 사용했다.

2.3.1 Zalesak’s disk

2차원 Zalesak’s disk는 Fig. 2와 같이 (1.0 m, 1.0 m)의 정사각형 도메인으로 구성되고, 초기에 Disk 물체가 (0.5 m, 0.75 m)에 위치해있다. Disk 물체는 Fig. 3과 같이 0.15 m의 반경을 가졌으며, 중심에서 높이 0.1 m, 두께 0.075 m의 직사각형 모양으로 뚫려있다.

Fig. 2. 
Computational domain for Zalesak’s disk.

Fig. 3. 
Schematic of Zalesak’s disk, unit m.

Fig. 4. 
Schematic of 3-D deformation, unit m.

계산 도메인에 전체에 적용된 2차원 강제 속도장은 다음과 같다.

u=-2πy-0.5υ=2πx-0.5

(14)

이 속도장으로 인해 Disk 물체는 그 모양을 유지하면서 반시계 방향으로 회전하고, T = 1.0초 뒤에 초기 위치로 돌아온다.

2.3.2 Zalesak’s disk 결과

Fig. 5는 격자 크기에 따른 Zalesak’s disk의 결과를 각각 t = 0, T/3, 2T/3, T에서의 F = 0.5의 선 컨투어를 나타낸 그림이다. 시간이 지남에 따라 LS와 VOF의 이송에 대한 오차가 쌓이기 때문에 표면의 변형이 심해진다. 또한 격자가 조밀해질수록 오차가 작아진다.

Fig. 5. 
Results for Zalesak’s disk along the time.

Fig. 6은 이론적인 Disk의 표면 분포(검정 선)와 t = T에서의 표면 분포(빨강선)를 확대 도시한 그림이다. 계산 결과를 이론적인 분포와 비교해 봤을 때 LS와 VOF 이송은 이동하는 방향 쪽으로 과하게 예측하고, 격자가 조밀해질수록 이론적 분포와 가까워진다. 또한 표면의 기울기가 급변하는 구간에서 변형이 심해진다.

Fig. 6. 
Results for Zalesak’s disk at F=0.5; Black line is the exact results, Red line is the calculation results.

Fig. 7은 격자에 따른 상대적 보존 오차의 그래프이다. 여기서 상대적 보존 오차는 다음과 같다.

Emass=∫Fxi,tdxi-∫Fxi,0dxi∫Fxi,0dxi

(15)

Fig. 7. 
The relative conservation errors for Zalesak’s disk.

상대적 보존 오차는 Eq. 8로 인해 시간에 따라서 진동하게 되고, 격자가 조밀해 질수록 진동 값이 작아진다. 최대 상대적 보존 오차는 h = 1/64에서는 0.0616, h = 1/128에서는 0.0427, 그리고 h = 1/160에서는 0.0316으로 나타났다.

2.3.4 액체의 3차원 변형 계산 결과

Fig. 8은 격자 크기에 따른 결과를 각각 t = 0, T/2 T에서의 F = 0.5의 iso-surface 컨투어를 나타낸 그림이다. 변형의 최대지점인 T/2 에서의 결과를 보면, 격자가 조밀할수록 나타낼 수 있는 표면의 크기가 작아지기 때문에 표면의 찢어짐이 완화되고 상대적 보존 오차가 작아지기 때문에 T에서의 표면이 초기 표면과 유사해 진다.

Fig. 8. 
Three dimensional deformation along the time.

Fig. 9(b)는 (a)에서 나타낸 위치에서의 VOF 컨투어 모습이고, Fig. 9(c)는 초기 표면(검정 선)과 T에서의 표면(빨간 선)을 비교한 그림이다. Fig. 9(b)를 보면 과소 오차에 의해 물체 내부에 VOF 값이 1보다 작아지는 곳이 발생하고, 이 오차의 크기는 격자가 조밀해 질수록 작아진다. Fig. 9(c)에서도 격자가 조밀해 질수록 에서의표면 모양이 초기 표면의 모양인 구형의 모양과 유사해짐을 알 수 있다.

Fig. 9. 
VOF contours shape for different grid sizes.

Fig. 10은 t = T에서 격자에 따른 VOF=0.01의 iso-surface 컨투어 모습이다. Fig. 10(a)를 보면 물체 외부에 VOF가 0보다 큰 곳이 생성되었고, 이는 물체가 변형함에 따라 표면을 잘 예측하지 못하여 과도 오차가 발생되었기 때문이다. 이 오차는 격자가 조밀해 질수록 작아지고, Fig. 10 (c)의 160³ 격자에서는 잘 나타나지 않는다.Fig. 11은 격자와 시간에 따른 상대적 보존 오차를 나타낸 그래프이다. Zalesak‘s disk의 경우와 마찬가지로 진동하는 모습을 보인다. 최대 상대적 보존 오차는 h = 1/64에서는 0.0604, h = 1/128에서는 0.0458, 그리고 h = 1/160에서는0.0232로 나타났다. Zalesak’s disk와 3차원 변형 문제에서의 최대 상대적 보존 오차는 Table 1에 나타나 있다.

Fig. 10. 
Three dimensional VOF iso surface for different grid sizes.

Fig. 11. 
The relative conservation errors for liquid 3-D deformation.

Table 1. 
The maximum relative conservation error of each grid size.

h	1/64	1/128	1/160
Case
Zalesak’s disk	0.0616	0.0427	0.0316
3D deform.	0.0604	0.0458	0.0232

비압축성 방정식이 결합된 2상 유동 문제를 검증하기 위해 Rayleigh-Taylor 불안정성을 계산 하였다. Fig. 12와 같이 (1.0 m, 4.0 m)의 계산 도메인을 사용하였으며, 128 × 512의 격자로 구성했다. 경계조건은 윗면과 아랫면에 No-slip의 벽면 조건을, 양 옆에는 Symmetric 조건으로 설정했다. 중간의 표면을 기준으로 위와 아래에 밀도가 각각 1.205 kg/m³, 0.1694 kg/m³인 유체를 위치했고, 점성계수는 두 유체 모두 0.0031 Ns/m²이다. 또한 유체 간 표면장력은 무시했다. 초기 속도 장은 0이고 고밀도 유체가 중력으로 인해 아래로 내려오게 된다. 초기 표면의 위치는 다음과 같다.

Fig. 12. 
The computational domain of Rayleigh-Taylor instability, unit m.

Fig. 13은 t = 0.2, 0.4, 0.6, 0.8 sec 에서의 밀도 컨투어 모습이다. 시간이 지날수록 밀도 차이로 인한 Rayleigh-Taylor 불안정성이 발생하는 모습을 관찰할 수 있다.

Fig. 13. 
Density contours of Rayleigh-Taylor instability.

Chandrasekhar[21]는 Rayleigh-Taylor 불안정성을 선형적으로 분석하여 다음과 같은 함수를 제안했다.

여기서 Λ는 시간 t에서의 진폭이고, Λ₀은 초기 진폭, 그리고 α̂는 성장률이다. Fig. 14는 시간에 따른 진폭의 분석 값과 수치결과 값을 비교한 그래프이다. Chandrasekhar의 분석에서 성장률 값은 5.05의 일정한 값으로 나타나지만 본 연구에서는 3.46에서 7.50까지 넓은 범위의 값으로 나타났다.

Fig. 14. 
Disturbance amplitude versus time.