하이브리드 입자-격자 방법에서의 압력장 계산

doi:10.3744/SNAK.2014.51.3.328

OA학술지
Journal of the Society of Naval Architects of Korea

하이브리드 입자-격자 방법에서의 압력장 계산 Computation of Pressure Fields for a Hybrid Particle-Mesh Method

DOI : 10.3744/SNAK.2014.51.3.328
Publish: Journal of the Society of Naval Architects of Korea Volume 51, Issue4, p328~333, 20 Aug 2014

ABSTRACT

A hybrid particle-mesh method based on the vorticity-velocity formulation for solving the incompressible Navier-Stokes equations is a combination of the Vortex-In-Cell(VIC) method for convection and the penalization method for diffusion. The key feature of the numerical methods is to determine velocity and vorticity fields around a solid body on a temporary grid, and then the time evolution of the flow is computed by tracing the convection of each vortex element using the Lagrangian approach. Assuming that the vorticity and velocity fields are to be computed in time domain analysis, pressure fields are estimated through a complete set of solutions at present time step. It is possible to obtain vorticity and velocity fields prior to any pressure calculation since the pressure term is eliminated in the vorticity-velocity formulation. Therefore, pressure field is explicitly treated by solving a suitable Poisson equation. In this paper, we propose a simple way to numerically implement the vorticity-velocity-pressure formulation including a penalty term. For validation of the proposed numerical scheme, we illustrate the early development of viscous flows around an impulsive started circular cylinder for Reynolds number of 9500.

KEYWORD

압력장 , 와도-속도-압력 정식화 , 보텍스인셀 방법 , 패널티화 방법

본문

Collapse all

1. 서 론

라그랑지안 보텍스 방법(Lagrangian vortex method)은 기본적으로 Navier-Stokes 방정식을 따라 와도 입자의 생성, 이송, 확산과정을 추적하는데, 세부 알고리즘에 따라 다양한 종류로 나누어진다. 보텍스 방법은 초기에 Lighthill (1963)과 Bachelor (1967)에 의해 연구가 되었으며, Sarpkaya (1989)은 많은 연구자들에 의해 개발된 다양한 보텍스 방법(Lagrangian or mixed Lagrangian-Eulerian schemes, the Biot-Savart law or the vortex-in-cell(VIC) methods)을 소개하였다. 와도를 기저로 하는 보텍스 방법은 비관성 및 관성 좌표계에서의 지배 방정식 형태가 같다. 또한 와도가 나타나는 물체 근처의 한정된 계산영역을 택할 수 있으며, 무한원방 조건이 자동으로 만족된다는 점에서 비정상 유동해석에 적합하다. 또한 보텍스 방법은 유동의 물리현상에 기초한 비교적 간단한 알고리즘으로 구성되기 때문에 복잡한 비정상 유동 해석에 유리하다. 이산화된 와도장과 라그랑지안 개념으로 표현된 와도 이송방정식을 기본으로 하기 때문에 보텍스 방법은 와도값을 가지는 영역만을 대상으로 계산이 이루어진다. 보텍스 방법은 압력항이 없는 와도와 속도만으로 표현된 지배방정식을 다룬다. 이는 압력항의 연성문제가 자연스럽게 해결되는 장점이 있지만 물체 주위의 압력을 구하기 위해서는 얻어진 와도장과 속도장을 가지고 압력장을 추가적으로 계산하여야 한다.

본 논문에서는 비정상 점성유동 해석을 위해서 VIC 방법과 penalization 방법이 결합된 하이브리드 입자-격자 방법을 사용한다. 보텍스 방법의 주된 문제는 속도-와도 관계식인 Biot-Savart 적분에 많은 시간이 소요되는 것과 수치기법에서 물체 표면에서의 표면 고착조건이 만족되도록 물체표면에서의 와도 값을 효과적으로 결정하는 것이다. 이를 해결하기 위해서 VIC 방법은 유동장에 위치한 와도 입자를 직교 좌표계로 보간(interpolation)한 후 정규격자상에서 고속 푸리에 변환(fast fourier transform, FFT)를 이용하여 와도의 이동 속도를 계산한다. 또한, Penalization 방법은 표면에서의 고착조건(no-slip condition)과 비침투 조건(no-penetration condition)을 동시에 만족시키면서 와도의 확산항 계산하는 방법으로 Angot (1999)가 penalty 항에 대한 수학적으로 증명한 이후 많은 수치기법에서 사용되고 있다. VIC 방법과 Penalization 방법은 병렬계산의 도입과 더불어 보텍스 방법의 계산 시간을 단축하는데 매우 효과적이다. 본 논문에서는 이러한 수치기법을 기반으로 와도장과 속도장을 구하며, 얻어진 결과로 부터 압력장을 계산하는 수치적 방법을 제안하고자 한다. 2장에서 지배방정식과 hybrid particle-mesh방법에서 사용되는 수치기법을 간단히 소개하며, 3장에서는 압력장 계산을 위한 수치방법을 설명한다. 제안된 수치기법의 검증을 위한 유동해석 결과는 4장에서 보여준다. 본 논문에서는 순간 출발하는 2차원 실린더에 대한 점성유동해석 수행하여 제안된 방법을 검증하고자 한다.

2. 지배 방정식

2차원 비압축성 점성 유동장에서 보텍스방법은 이산화된 와도장을 표현하는 보텍스 입자를 기본 요소로 한다. 와도장은 N개의 보텍스 입자(순환세기 Γ=ωh²)와 순정함수(𝜂_𝜖)를 이용하여 식 (1) 과 같이 표현된다.

보텍스 입자의 성장과 이동은 식 (2)와 (3)에서와 같이 Lagrangian 개념으로 표현된다.

2.1 Vortex-In-Cell 방법

점성유동자의 속도(u=U_∞+u_ω)는 Helmholtz 분해로 표현하면 유입 속도 U_∞ 와 와도가 유기하는 회전 속도 u_ω로 구성된다. 회전속도 성분을 얻기 위해서 유량함수와 와도 간의 Poisson 방정식 ∇²𝜓=−ω을 먼저 계산하여 유량함수(𝜓)를 얻은 후, 관계식 에 의해 속도를 구한다. 또한, Dirichlet 경계조건을 사용하기 위해서 경계에서의 유량함수는 Green 함수의 해를 이용하여 구하였다. VIC 방법에서 유량함수와 속도의 계산은 정규격자계에서 수행되며, 와도 값은 식 (4)에 의해 정규격자계로 보간된다.

여기서 아래첨자 g와 p는 각각 와도 격자와 입자 상의 물리량을 의미하며, W 는 보간 함수(interpolation kernel function)으로 kernel을 사용하였다.

정규격자상에서 구해진 속도 성분은 다시 식 (4)를 이용하여 입자 위치에서의 속도를 구한다.

2.2 Penalization 방법

표면 고착조건을 만족시키기 위한 penalty 항이 추가된 Navier-Stokes 방정식은 식 (5)와 같다.

여기서 u_s 는 물체의 속도이다. 𝜒 는 mask 함수로서 물체 안에서는 1이고 물체 밖에서는 0이 된다. λ는 penalty 상수로서 매우 큰 값이어야 하며, Explicit 방법에서 λ는 1/Δt로 정의한다. 하지만 λ가 매우 큰 값을 갖기 위해서는 Δt가 과도하게 작아야하기 때문에 implicit penalty 항을 이용하여 와도 이송방정식을 다음과 같이 이산화한다.

확산항은 보텍스입자의 성장을 표현하며, Penalty 항은 물체표면에서의 고착조건을 만족하도록 물체 근처의 입자세기를 변경한다. 보텍스 입자의 속도 계산과 같이, 입자세기는 모두 정규격자계에서 구해진다. 변경된 입자세기는 입자위치로 다시 보간된다. 보텍스입자들은 최종적으로 얻은 속도와 세기를 갖고 이동하게 되며, 이때 시간전진은 midpoint predictor-corrector 방법을 사용하였고, Penalization parameter는 λ=10⁸을 사용하였다.

3. 압력 방정식

2차원 비정상, 비압축성 점성 유동장에서 Poisson 방정식 형태의 압력방정식은 penalty 항이 포함된 Navier-Stokes 방정식으로부터 얻어지며 식 (7)과 같다.

여기서 H 는 전체 압력수두를 나타내며 식 (8)과 같이 정의할 수 있다.

p_∞ 와 U_∞ 는 무한 원방에서의 압력과 속도를 의미한다. 따라서, 무한원방에서의 전체 압력수두는 0이 된다. 전체 압력수두는 VIC 방법에서와 같이 고속 푸리에 변환(fast fourier transform, FFT)을 기반으로 한 Poisson solver를 이용하여 계산할 수 있다. 식 (7)의 penalty항에서 penalization 상수 λ의 결정이 중요하게 된다. 와도 이송방정식에서 λ는 Δt 과도하게 작아 지는 것을 막기 위해서 implicit 방법을 사용하며, 이때의 λ는 Δt상관없이 매우 큰 값(예, 10⁸)을 사용한다. 하지만 압력방정식에서 같은 값의 λ를 사용할 수 없다.

문제를 단순화 하기 위해 물체속도가 없는 경우를 생각한다면, 식 (7)은 식 (9)과 같이 표현된다.

여기서 λ′는 와도 이송 방정식에서의 λ와의 혼동을 피하기 위해서 사용된 penalization 상수이다. 이산화된 와도 이송방정식에서 penalization 상수가 무한히 큰 수라고 가정하면, λ𝜒(1+λΔt𝜒) ≈ 1/Δt가 된다. 이에 따라 λ′는 1/Δt로 생각할 수 있다. 보텍스 방법에서는 안정성 조건(stability condition)으로 classical CFL조건(Δt＜CΔx)이 아닌 νΔt/h²=0(1) 과 |ω|h²/ν = 0(1) 을 사용하며, 두 조건에서 |ω|Δt = 0(1)를 유도할 수 있다. 따라서 λ′ ≈ 1/Δt ≈ |ω|이며, 결국, 식 (9) 또는 식 (7)의 오른쪽 두 항의 물리량의 차수(order)가 유사해 짐을 알 수 있다. 본 논문에서 가정된 penalization 상수 λ′ = 1/Δt는 4장에서 수치적으로 검증한다.

4. 유동 해석

본 연구에서는 하이브리드 입자-격자 방법을 이용하여 순간적으로 출발하는 2차원 실린더 주위의 점성유동을 해석해 보았다. 이 문제는 단순하고 뭉뚝한 형상이면서도 복잡한 유동 특성을 포함하고 있기 때문에 수치기법 검증을 위한 기준 문제로서 적합하다. 레이놀즈 수 9500에 대한 유동해석을 수행하였으며, 유동 수치해석 조건은 Table1과 같다.

[Table 1] Computational parameters for Re=9500

Computational parameters for Re=9500

Fig. 1은 순간적으로 출발하는 2차원 실린더 주위의 와도장을 나타낸 것이다. 계산된 와도장은 Cottet, et al. (2010)의 결과와 잘 일치 함을 확인하였다. Fig. 2는 계산된 와도장을 바탕으로 계산된 속도장을 보여주며, T=2와 4에서의 속도장은 Loc and Bouard (1985)의 결과와 비교되었고, Bouard and Coutanceau (1980)가 실험에서 관찰한 α와 β 현상이 수치적으로 잘 모사되었음을 확인 할 수 있다. Bouard and Coutanceau (1980)는 α와 β 현상을 Fig. 3과 같이 정의하고 있다. α현상은 2차 와류(secondary vortex)가 두 개로 나누어지는 것을 말하며 레이놀즈 수 500 이상의 실린더 주위 유동에서 나타난다. 또한 β 현상은 초기 발생한 1차 와류(primary vortex)가 점차 사라지면서 새로운 기본적인 보텍스가 형성되는 것을 말한다. β 현상은 레이놀즈 수 3000 이상의 유동에서 관찰된다.

[Fig.1] Instantaneous vorticity distribution for Re = 9500; (a) T=1, (b) T=2, and (c) T=3. Note that the figures on the right-hand side column adapted from Cottet et al. (2010)

[Fig. 2] Instantaneous streamlines for an impulsively started circular cylinder: (a) T = 2 and (b) T = 4 compared with flow visualization results of Loc and Bouard (1985) for Re = 9500

[Fig. 3] Schematic representation of the secondary eddy pair: (a) phenomenon α (left); (b) phenomenon β (right). Note that the figures are reprinted from Bouard and Coutanceau (1980)

잘 모사된 와도장과 속도장을 바탕으로 2차원 원형 실린더 주위의 압력장을 구하였으며, 계산된 압력장은 Fig. 4와 같다. Penalization 방법의 특성상 물체 내부의 압력장도 같이 구해지며, 물체 경계를 기준으로 물체 내부와 외부의 압력변화가 부드럽게 이루어지고 있음을 볼 수 있다. 하지만 물체내부의 압력장은 실제 존재하는 것은 아니며 계산 상에서만 존재하는 압력 분포이다. Fig. 5는 물체표면에서의 압력분포를 보여준다. T=1, 2, 3에서의 물체 표면 압력계수는 유한체적법(finite volume method)으로 해석한 Suh and Kim (1999)의 결과와 매우 잘 일치하고 있다.

[Fig. 4] Pressure coefficient distribution for Re = 9500

[Fig. 5] Surface pressure distribution at several instants for Re = 9500 compared with numerical results of Suh and Kim (1999)

계산영역이 충분히 크다고 가정하기 때문에 식 (9)에서 경계조건 H = 0을 사용하였다. 하지만 균일한 경계조건(homogeneous boundary condition)에 의한 해는 계산영역에 따라 다르게 나타날 수 있다. 이러한 단점을 개선하고자 적분식을 이용한 비균일 경계조건(non-homogeneous boundary condition)을 제안하며, 컴퓨터 메모리의 한계로 인하여 계산 영역에 대한 검증은 레이놀즈수 550의 유동에 대해서 수행하였다.

계산 영역의 경계에서의 Navier-Stokes 방정식은 penalty 항이 사라지게 되므로 식 (10)과 같이 표현된다.

계산영역은 모든 보텍스 입자를 포함하도록 선택하기 때문에 경계에서의 와도값은 0이 되므로 식 (10)에서 좌변의 두번째 항은 사라지게 된다. 식 (10)의 양변에 적분을 취하면 식 (11)을 얻을 수 있으며, 식 (12)와 같이 이산화된다.

여기서 n = 1, 2, 3, …, N 이고, N 은 경계에서의 격자수이다. dl은 격자점에서의 접선벡터을 의미한다. 로 가정하면, 전체 압력수두 값은 계산영역의 경계를 따라서 식 (12)에 의해서 계산될 수 있다. Fig. 6은 균일 경계조건과 적분방법에 의한 비균일 경계조건에 의해 구해진 압력장을 보여준다. 적분방법에 의해 계산된 압력장이 균일한 경계조건에 의해 구해진 것보다 계산영역의 영향이 작음을 알 수 있다. 하지만, Fig. 7에서와 같이 물체 경계에서의 압력계수 값의 차이가 매우 작다. 계산영역 [1024 × 1024]에 대해 균일한 경계조건을 사용한 결과를 제외한 계산결과에서 표면 압력계수의 차이는 0.2％이내 이며, 균일한 경계조건의 계산영역 [1024 × 1024]의 결과는 다른 결과와 비교하였을 때 표면 압력계수의 차이가 3％이내 이다. 적분방법에 의한 경계값의 계산에서 적분 시작점의 압력을 무한원방의 압력과 같다고 가정하고 있으나, 적분방법에 의한 비균일 경계조건은 균일 경계조건에 비해 상대적으로 계산 영역에 대한 영향이 작기 때문에 계산영역이 과도하게 커지는 것을 막을 수 있다.

[Fig. 6] Pressure distribution for Re = 550 at T=3

[Fig. 7] Surface pressure coefficients at different domain sizes

5. 결 론

본 연구에서는 와도를 기본 변수로 하는 하이브리드 입자-격자 방법에서 압력장을 계산하기 위한 수치기법을 제안하였다. 제안된 방법은 순간적으로 출발하는 원형실린더에 대한 2차원 유동해석 결과를 통하여 검증하였다. 또한 논문에서 제시한 적분방법에 의한 경계값의 추정은 상대적으로 작은 계산영역을 취할 수 있기 때문에 압력장 계산 수행에 있어서 효율적이다. 와도-속도-압력 정식화를 기반으로 한 하이브리드 입자-격자 방법은 비교적 간단한 알고리즘으로 구성되며 물체와는 상관없는 정규 직교 격자계를 사용하기 때문에, 해석코드 개발이 비교적 쉽고 복잡한 형상을 갖는 물체에 대한 적용이 유리할 것으로 판단된다.

참고문헌

1. Angot P., Brunear C.H., Fabrie P. 1999 A penalization method to take into account obstacles in incompressible viscous flows [Numerische Mathematik] Vol.81 P.497-520
2. Batchelor G.K 1967 An introduction to fluid dynamics
3. Bouard R., Coutanceau M. 1980 The early stage of development of the wake behind an impulsively started circular cylinder for 40 ＜ Re ＜ 104 [Journal of Fluid Mechanics] Vol.101 P.586-607
4. Cottet G.H., Gallizio F., Magni A., Mortazavi I 2010 A vortex immersed boundary method for bluff body flows [Proceedings of 3rd Joint US-European ASME Fluids Enigineering Summer Meeting, Sympoisum on Development and Applications of Immersed Boundary Methods]
5. Lighthill M.J., Rosenhead J. 1963 Introduction, boundary layer theory, Laminar boundary layers
6. Loc T.P., Bouard R. 1985 Numerical solution of the early stage of the unsteady viscous flow around a circular cylinder: a comparison with experimental visualization and measurements [Journal of Fluid Mechanics] Vol.160 P.93-117
7. Sarpkaya T 1989 Computational methods with vortices: the 1988 Freeman scholar Lecture [Journal of Fluids Engineering] Vol.111 P.5-52
8. Suh J.C., Kim K.S. 1999 A vorticity-velocity formulation for solving the two-dimensional Navier Stokes equations [Fluid Dynamics Research] Vol.25 P.195-216

이미지 / 테이블

[ ]
[ ]
[ ]
[ ]
[ ]
[ ]
[ ]
[ ]
[ ]
[ ]
[ Table 1 ] Computational parameters for Re=9500
[ Fig.1 ] Instantaneous vorticity distribution for Re = 9500; (a) T=1, (b) T=2, and (c) T=3. Note that the figures on the right-hand side column adapted from Cottet et al. (2010)
[ Fig. 2 ] Instantaneous streamlines for an impulsively started circular cylinder: (a) T = 2 and (b) T = 4 compared with flow visualization results of Loc and Bouard (1985) for Re = 9500
[ Fig. 3 ] Schematic representation of the secondary eddy pair: (a) phenomenon α (left); (b) phenomenon β (right). Note that the figures are reprinted from Bouard and Coutanceau (1980)
[ Fig. 4 ] Pressure coefficient distribution for Re = 9500
[ Fig. 5 ] Surface pressure distribution at several instants for Re = 9500 compared with numerical results of Suh and Kim (1999)
[ ]
[ ]
[ ]
[ Fig. 6 ] Pressure distribution for Re = 550 at T=3
[ Fig. 7 ] Surface pressure coefficients at different domain sizes