Flow Analysis around a Wing Section by a Piecewise Linear Panel Method

부분선형 패널법을 이용한 2차원 날개단면 주위 유동 해석

Abstract

Panel methods are useful tools for analyzing fluid-flow around a wing section. It has the advantage of fast and accurate calculation, compared to other CFD Methods such as RANS solvers. This paper suggests a piecewise linear panel method in order to improve accuracy of existing panel methods by changing the piecewise constant singularity strength to linear singularity strength(for dipole strength). The piecewise linear panel method adopts the linear distribution of singularity strength, while control point is located at the node of each panel. Formulation of the piecewise linear panel method is given, and some calculation results are shown for typical wing sections.

1. 서 론

프로펠러 주위의 비점성 유동 해석을 위해 전통적으로 패널법을 널리 사용해오고 있다. 이는 패널법을 이용한 유동해석의 계산 시간이 매우 짧으며 높은 정도를 보이는 장점 때문이라 할 수 있다. 하지만 최근 컴퓨터 성능의 발전에 따른 전산유동해석기법(CFD, Computational Fluid Dynamics)이 점차 다양한 분야에서 많이 활용되고 있으며 그 역할이 커질 전망이다. 그러나 전산유동해석기법은 격자의 질과 종류에 따라 해석에 많은 영향을 끼쳐 그 정도를 확인하는데 제한이 있으며, 해석 시간이 패널법과 비교하여 상당히 오래 걸린다는 단점이 있다. 그리하여 패널법은 여전히 유용하고 효과적인 수단으로서의 가치가 있다 (Park, 2015).

본 연구에서는 2차원 날개단면 주위 유동해석에 대해 기존의 패널법을 부분선형 패널법으로 보완 및 발전시키고자 하였다. 즉 기존 패널법에서는 특이점의 세기가 패널에 걸쳐 일정하다고 가정하였다. 하지만 본 연구에서는 Jack Moran (1983)이 제안한 방법을 기반으로 패널에 걸쳐 선형적으로 변하는 다이폴 특이점의 세기분포를 가진 부분선형 패널법을 개발하였다. 쏘스(source)에 대한 부분은 기존과 같이 일정한 세기를 갖는다.

부분선형 패널법에서는 특이점의 세기가 노드(Node)에서 연속되므로 기존 방법에서 특이점 세기분포의 단차로 인한 오차를 줄이고 2차원 날개단면 주위 유동해석의 정도 향상을 기대하였다. 또한, Moran의 방법을 수정하여 계산의 편리성을 도모하고자 하였다.

 

2. 패널법의 지배방정식

2.1 가정과 경계조건

패널법 문제에서는 유동장 V, 그 경계면을 S, 유동장 방향을 향하는 법선벡터 을 정의한다. 경계면 S는 물체의 표면(body surface) SB , 후류면(wake surface) SW, 그리고 물체표면과 후류면 주변의 외부면(outer control surface) S∞으로 구성된다. 또한 물체를 향하는 방향으로 유입속도 가 있다.

유동장내의 유체는 비압축(incompressible), 비점성(inviscid) 그리고 비회전성(irrotational)으로 가정하며, 따라서 라플라스 방정식(Laplace equation)을 만족하는 섭동속도 포텐셜(perturbed velocity potential)이 존재한다고 가정한다.

경계치 문제에서는 경계면 S에서의 경계조건을 다음과 같이 나누어 설정 할 수 있다.

물체표면 SB 에서 운동학적 경계조건을 만족해야한다

후류면 SW는 두께가 없다고 가정하며, 운동학적 경계조건으로 후류상의 상,하면의 법선속도는 같고, 동역학적 경계조건으로 후류면에 작용하는 힘은 없다. 또, 운동학적 경계조건으로부터 법선속도가 같다면, 즉 포텐셜의 법선 방향 변화량이 0 이며, ∆ϕ= const로 상하면의 포텐셜 점프(차이)가 있다는 것을 내포하고 있다.

Fig. 1Notation for a general body for the application of Green's Theorem

• 후류면을 가로지르는 포텐셜 점프는 물체 주변의 순환(circulation)의 크기와 같으며 후류면을 따라 일정하다.

• 날개 뒷날에서의 속도는 유한하다는 Kutta 조건을 만족한다.

• 물체의 외부면 S∞에서는 물체에 의한 섭동 속도는 물체로부터 무한대의 거리에 있는 곳에서 사라진다.

2.2 기본 방정식

우선 패널법에 쓰이는 형상을 정의하면 Fig. 2와 같다. 날개 형상을 N개의 패널로 나누면 N+1개의 노드가 생성된다. 노드와 패널의 순서는 날개끝단 아래 부분부터 시작하여 형상의 앞날을 돌아 다시 끝단의 윗면 순으로 정의된다.

Fig. 2Representation of a hydrofoil with straight line segments

그리고 각 패널의 국부좌표는 Fig. 3과 같이 정의하였다.

Fig. 3Nomenclature for local coordinate systems

ith패널은 ith노드와 i+1th노드 사이에 놓인다. 그리고 패널의 x축에 대한 각도를 θ로 둔다. 이러한 가정 하에 θ에 대한 sin, cos값은 식 (8)과 같다.

그리고 법선과 접선벡터는 식 (9)와 같이 정의하였다

제어점(control point)은 Fig. 4에서 보듯이 각 패널의 절점(node)에 위치한다.

Fig. 4Local panel nomenclature

이러한 좌표계의 정의아래, 유동장내 물체표면에 위치한 제어점 p에서 특이점 q에 의한 포텐셜은 다음과 같이 물체 표면에 위치한 다이폴, 쏘스, 그리고 웨이크에 의해 결정된다.

where,

ϕ(p) = Perturbation potential in V

p(x,y) = Field point where induced potential is calculated

q(ξ,ƞ) = Singularity point

R(p;q) = Distance between point p and q

= Normal derivative with respect to the point q

ϕ(q) = Dipole Strength at q

= Source Strength at q

이때, 쏘스 세기, 는 경계조건 식 (2), 즉 으로부터 알 수 있다.

따라서 식(10)은 물체표면에서 다이폴 세기, ϕ(q)를 미지수로 하는 적분방정식문제라 할 수 있다. 그리고 식(10)의 이산화는 미지수 ϕ(q)에 대한 선형 방정식이 된다.

 

3. 부분 선형 패널법

3.1 소개

기존의 패널법에서는 Fig. 5의 (a)처럼 적분 방정식의 이산화식을 만족하는 제어점(Control Point)은 패널의 중앙위치하며, 패널 내부 특이점 강도가 패널에 걸쳐 일정하게 분포한다고 가정하였다.

Fig. 5Influence of constant and linear strength panels

본 연구의 부분선형(piecewise linear) 패널법에서는 단일 패널에 걸쳐 특이점의 강도가 선형적으로 변함을 의미한다. 이는 Fig. 5(b)처럼 절점에서 특이점의 세기는 인접한 두 패널의 선형적으로 변하는 특이점의 강도에 의한 영향으로 결정된다고 가정하였다.

기존의 방법은 패널마다 특이점이 일정한 크기를 가지므로 인접한 패널과의 특이점세기의 단차가 필연적으로 발생한다. 하지만 부분선형패널법에서는 특이점의 세기가 인접한 두 패널에 의해 결정되므로 특이점분포가 단차 없이 연결됨을 알 수 있다.

이러한 방법으로 인한 기대효과는 기존의 패널법이 가지는 특이점 세기에 대한 오차를 줄 일 수 있으리라 예상되었다.

3.2 특이점의 이산화

식 (10)은 이산화에 의하여 다음과 같은 선형 대수방정식으로 치환된다.

여기서,

식 (11)과 같은 기존의 균일강도 패널법에 선형강도를 가지는 다이폴을 적용하게 되면 다음과 같다.

여기서,

Aij = DRij−1 + DLij if i = 1,2,3,,,, N j = 2,3,4,,, N

Aij = DLij i f j = 1 = DRij if j = N + 1

Aij = Aij − Wi if j = 1 = Aij + Wi if j = N + 1

위의 수식에서 DL과 DR은 다음 그림과 같은 강도를 가지는 다이폴 영향계수를 의미하고, 그 계산은 appendix에 보였다.

Fig. 6DR dipole influence coefficient

Fig. 7DL dipole influence coefficient

균일강도 패널법에서는, 제어점과 특이점의 위치가 같을 때 (즉, i = j) 다이폴에 의한 Self-induction이 발생하고 이는 다이폴의 영향 계수에 의해 그 크기는 π가 된다.

반면, 선형강도 패널법에서는 제어점의 위치가 절점에 위치함으로 각 절점에서의 self-induction을 고려하여야 한다. 따라서 본 논문에서는 Jack Moran이 제안한 방법을 사용하였다. Jack Moran은 각 절점에서의 self-induction은 절점에 인접한 두 패널의 각도로 표현 할 수 있다고 보았다.

Fig. 8에서 제어점 P가 ith 절점이 되면 self-induction은 다음과 같이 표현할 수 있다.

Fig. 8Self-induction for linear dipole strength panel method

위와 같이 계산된 self-induction을 식 (12)에 대입하면, i = j (i, j = 2, 3 ..., N) 일 때 Aij = θself가 된다. 반면 Ai,j = A1, 1,그리고 Ai,j = A1 , N + 1일 경우, 각 꼬리 날개 각도의 반이 된다.

3.3 Kutta 조건

부분선형패널법에서는 N개의 방정식과 N+1개의 미지수가 있기 때문에 방정식을 풀기위한 1개의 식이 추가적으로 요구된다.

Fig. 9와 같이 Kutta 조건을 세워 적용하였다.

Fig. 9Self-induction of Node 1 and N+1

1번 패널제어점에 유기되는 접선방향 속도(tangential velocity)와 N번 패널제어점에 유기되는 접선방향 속도의 크기가 같도록 하는 Kutta 조건을 세웠다.

Fig. 10Notation for Kutta Condition at the trailing edge

여기서,

 

4. 계산 결과 및 결론

4.1 계산 결과

기존의 Jack Moran이 제안한 선형 패널법에선 N개의 방정식과 N개의 미지수를 두고 풀었던 것과는 달리 본 논문에서는 Kutta 조건을 적용하여 N+1개의 방정식과 N+1개의 미지수를 이산화 방정식을 풀어서 구하였다.

1) 수식의 검증을 위해 원 단면에 대하여 입사각이 0°일 경우와 90°일 경우에 대해 해석해에 대해 패널 개수가 10, 20 ,40 그리고 80개로 증가함에 따른 포텐셜의 수렴도를 확인 하였다.

Fig. 11과 12에서 보듯이 패널의 개수가 증가함에 따라 수렴하고 있음을 알 수 있고, 이를 통해 수식을 신뢰할 수 있다고 보인다.

Fig. 11Convergence test of potential at 0 degree

Fig. 12Convergence test of potential at 90 degree

2) 수식의 정밀도를 보이기 위해 Karman-Trefft z 단면에 대해 해석해, Moran에 의해 제안된 방법 그리고 본 논문에서 제안한 방법을 패널의 개수가 80개 일 때의 결과 값들을 비교하였다.

Fig. 13에서 보듯이 Moran이 제안한 방법과 본 논문에서 제안한 방법의 정밀도가 같음을 알 수 있다. 다만 두 방법 모두 해석해와는 차이를 보이고 있으며, 이는 날개의 꼬리부분이 매우 얇아서 생기는 오차에 의한 것이라고 추측된다.

Fig. 13Comparison of potential values for analytic, Moran and present methods for a K-T section (No.of panel = 80, tau = 1.95, x = 0.1, y = 0.1, AOA = 0deg)

3) 선형 강도 패널법의 정밀도를 보이기 위해 NACA0012에 대해 입사각이 90°일 때 −Cp 값을 균일 강도 패널법과 비교하였다.

Fig. 14에서 보듯이 균일강도 패널법에 비해 선형강도 패널법이 날개의 앞에서 피크 값을 더 잘 표현하고 있음을 알 수 있다.

Fig. 14Comparison of −Cp values for constant and linear panel methods for NACA0012 section. (No.of panel = 80, AOA = 90deg)

4.2 결 론

균일강도 패널법은 패널에 걸쳐 강도가 일정하다고 보았다. 그로인해 각 패널 간 강도의 불연속이 발생하게 된다. 본 논문에서는 균일강도 패널법에서 발생하는 불연속을 선형강도 패널법을 이용해 해결하고자 하였다. 또한, Moran이 제안한 선형강도 패널법을 수정하여, N + 1th 절점의 계산 값을 이산화 방정식을 풀어 결과를 얻고자 하였다.

균일강도 패널법의 패널 간 강도의 불연속을 패널에 강도를 선형적으로 분포시켜 해결하여 계산결과의 정도를 향상시킬 수 있었다. 또한, Moran의 방법을 수정하여 계산결과의 정도는 유지하면서 N + 1th 절점의 값을 보간법을 사용하지 않고 이산화 방정식을 풀어서 계산할 수 있게 되었다. 추후 본 논문에서 제안한 방법을 3차원 프로펠러 형상 해석에 적용하면, 계산 결과의 정도가 향상 될 것이라고 예상된다.

 

Appendix

영향계수의 계산

A.1 좌표계

각 패널에 해당하는 국부 좌표계 (ξ,η)와 전체적인 전 좌표계 (x,y)에 대해 정의하였다.

Fig. 15The local and global coordinate system

xc , yc : Field point where induced potential is calculated tx , ty : Distance from the field point to the first panel point in local coordinate R1 : Distance from the field point to the first panel point R2 : Distance from the field point to the second panel point dx1, dy1 : Distance from the field point to the first panel point in global coordinate dx2, dy2 : Distance from the field point to the second panel point in global coordinate

법선, 접선 벡터를 정의하였다.

국부 좌표계에서 적절한 법선 벡터를 이용하여 전 좌표계로 바꿔줄 수 있다.

여기서,

Us : X-direction Velocity in Global Coordinate

Vs : Y-direction Velocity in Global Coordinate

Uk : ξ-direction Velocity in Local Coordinate

Vk : η-direction Velocity in Local Coordinate

A.2 영향 계수

A.2.1 균일 강도

영향 계수 계산을 위한 수식에 쓰인 값들을 정리하면 다음과 같다.

우선 특이점 쏘스에 의한 포텐셜은 다음과 같다.

다음으로 특이점 다이폴에 의한 포텐셜은 다음과 같다.

A.2.2 선형 강도

부분선형 패널법에서의 쏘스는 균일강도 패널법에서 사용한 영향 계수와 같으므로 다이폴에 의한 영향 계수의 계산만 다루도록 한다.

선형 강도는 패널에 걸쳐 시작점에서 끝점으로 갈수록 강도가 증가하는 것과 감소하는 것으로 구분할 수 있다. 우선 Fig. 6에서 보듯이 다이폴의 강도가 증가하는 영향 계수의 계산은 다음과 같다.

다음으로 Fig. 7에 보인 강도가 감소하는 영향 계수의 계산은 다음과 같다.