3.2. 특이점 분포의 이산화
기존의 저차 패널법에서 이산화식을 만족하는 제어점은 패널의 중앙에 위치하고, 패널에 걸쳐 다이폴의 강도가 일정하게 분포되어 있다고 가정하였다. 본 연구에서는 인접한 패널들간에 다이폴 강도의 단차를 제거하고자 패널에 걸쳐 다이폴의 강도가 선형적으로 분포한다고 가정하고, 제어점과 특이점은 각 패널의 꼭짓점에 위치한다. 이러한 방법은 2차원 단면에 대해 연구가 진행되어 검증된 바 있다 (Park et al., 2015).
경계치 문제의 해는 적분방정식 식 (10)을 만족하는 특이점의 강도를 결정하는 것으로 구성된다. 쏘오스의 강도는 운동학적 경계조건에 의해 결정되고 오직 다이폴의 강도를 구하는 문제로 귀결된다 (Kim, 2003).
적분방정식 식 (10)을 이산화 시키면 다음과 같이 표현된다.
∑j=1NnodeDijϕj+∑m=1MrWim(∆ϕ)m=∑j=1NnodeSij(∂ϕ∂n)j, i=1,2,...,Nnode | (11)
|
여기서, Nnode은 총 절점의 개수이고 Nnode=(Nc2+1)×Mr+1로 표현된다. ∂ϕ∂n은 운동학적 경계조건에 의한 쏘오스의 강도이고, (∆ϕ)m은 m번째 반경의 날개 뒷날에서의 포텐셜 점프이다. Dij는 날개에 분포하는 다이폴에 의한 영향계수이고, Sij는 날개에 분포하는 쏘오스에 의한 영향계수이다. Wim은 후류면에 존재하는 다이폴에 의한 영향계수이다.
이산화 식 (11)을 선형강도를 갖는 다이폴을 적용하여 선형대수 방정식으로 표현 하면 다음과 같다.
∑j=1NnodeAijϕj =∑j=1NnodeBij(∂ϕ∂n)j, i=1,2,…,Nnode | (12)
|
여기서,
Aij=Dij ifi≠j = Solidangle ifi=j
Aij=Aij-Wi if i=1 =Aij+Wi if i=Nc+1
여기서, Aij는 다이폴에 의한 영향계수이고, Dij는 날개에 분포하는 다이폴이 선형적으로 분포할 때 다이폴에 의한 영향계수이다. 입체각(solidangle)은 날개상에 분포하는 다이폴에 의한 자기유도 영향계수로 꼭짓점을 공유하는 패널들이 이루는 3차원 각으로 표현되고, 3.3절에 자세히 나타내었다. Wi는 후류면 상에 존재하는 다이폴에 의한 영향계수이다. Bij는 쏘오스에 의한 영향계수이다.
선형 강도의 분포는 Fig. 3에 보는 것처럼 패널의 각 꼭짓점에 단위강도의 다이폴이 존재하고, 그 외의 꼭짓점에서의 강도는 영이 된다고 가정하여 영향계수를 계산하였다.
Fig. 3
Linear dipole strength distribution on a triangular panel
Fig. 3에서 보인 선형 강도의 분포는 다음의 식으로 나타낼 수 있다.
μi(xi, yi, zi)=Axi+Byi+Czi+D, i=1,2,3 | (13)
|
여기서, μi는 i번째 꼭짓점에서의 다이폴 강도이고, A, B, C 그리고 D는 다이폴의 강도로부터 결정되는 계수이다. 패널 각 변의 방향벡터를 eli→로, 꼭지점의 위치를 vi(xi, yi, zi)로 각각 정의한다. 선형강도 분포에 상응하는 패널의 기하학으로부터 A, B, C를 계산할 수 있다.
∇μ∙eli=μi+1-μili=A(xi+1-xi)+B(yi+1-yi)+C(zi+1-zi)li | (14)
|
여기서, li는 패널 각 변의 길이이고, 패널의 꼭짓점과 각 변은 반시계 방향으로 정의 된다. 식 (14)를 Cramer’s rule을 이용하여 계수 A, B, C를 계산한다.
A=det(μi-μi+1,yi-yi+1,zi-zi+1)∆B=det(xi-xi+1,μi-μi+1,zi-zi+1)∆C=det(xi-xi+1,yi-yi+1,μi-μi+1,)∆ | (15)
|
계산된 계수 A, B, C를 이용하여 D를 구한다.
다이폴에 의한 영향계수를 계산하는 경우, 제어점의 위치와 특이점의 위치가 같으면 자기 자신에 의해 유도되는 영향계수로 정의된다. 기존의 저차 패널법에서는 제어점의 위치가 패널의 가운데에 존재하여 그 값이 2π로 일정하다. 그러나 당 논문과 같이 삼각형 패널을 채택하는 경우 제어점의 위치가 패널의 꼭짓점에 위치하게 되며, 자기유도 영향계수는 꼭짓점을 공유하는 패널들사이의 각으로 표현되며 이는 각 꼭짓점을 잇는 내각으로 구성되는 입체각(solid angle)으로 표현된다 (Lee & Shu, 1995).
후류면에서의 다이폴은 사각 패널상에 쌍선형(bilinear)적으로 분포하고 있다고 가정하고, 그에 의한 영향계수는 Suh 등(1992)이 제안한 방법에 의해 계산한다 (Fig. 4).
Fig. 4
Bilinear dipole strength distribution on a wake panel
패널은 ζ=0의 평면에 있고, ζ축의 방향이 패널법의 법선 벡터의 방향과 같도록 ξ-η 국소 직교 좌표계(local coordinate system)를 정의한다. ξ축과 η축 방향의 단위 벡터는 각각 e→ξ와 e→η로 나타낸다.
후류 패널상 다이폴의 쌍선형 강도 분포는 다음과 같이 나타낼 수 있다.
μi(ξi,ηi) = Aξi+Βi+Cξiηi+D, i=1,2,3,4 | (16)
|
여기서, μi는 다이폴의 강도이고, A, B, C 그리고 D는 꼭짓점에서의 다이폴 강도에 의해 구해진 상수이다. 후류면에서의 다이폴 강도는 후류면의 양변(S1, S2)에서 정의된다.
본 논문에서는 제어점이 꼭짓점에 있다고 가정하였다. 하지만 날개 뒷날에서 제어점의 위치가 패널의 꼭짓점에 존재하게 되면 kutta 조건을 적용함에 있어서 어려움이 존재하게 된다. 이를 회피하기 위해 Fig. 5와 같이 날개 뒷날에서의 제어점을 날개 앞날 방향으로 이동시켜 계산을 진행 하였다.
Fig. 5
Relocation of the trailing edge control points
제어점의 적절한 위치를 파악하기 위해 입사각이 5.73°일 때 단면 NACA0005를 갖는 원형날개에 대해, 제어점의 위치를 날개 뒷날 패널 길이의 10~50% 까지 10% 단위로 앞날 방향으로 이동시켜 압력 Kutta 조건을 만족시키기 위해 수행된 반복계산 결과를 비교 하였다. 결과의 비교는 반복계산의 횟수와 날개 뒷날에서의 반경방향 최대 압력계수차를 기준으로 하였으며, 결과를 Table 1에 보였다. Table 1에서 보면 반복계산 횟수는 끝점에서 멀어질수록 줄어들고 있음을 확인 할 수 있다. 압력계수차의 경우 멀어질수록 커지는 현상을 보이지만 그 차이는 크지 않고 50%일 경우에도 충분히 수렴했다고 판단된다. 따라서, 본 연구에서는 제어점을 날개 앞날의 방향으로 50% 이동시켜 계산을 진행 하였다.
Table 1.
Location of the control point
Location of control point
|
Number of iteration
|
Maximum delta CP |
0.1
|
104
|
0.00001
|
0.2
|
63
|
0.00009
|
0.3
|
42
|
0.00004
|
0.4
|
23
|
0.00041
|
0.5
|
23
|
0.00014
|
3.3. 입체각
본 연구에서는 제어점과 특이점의 위치가 동일한 경우의 다이폴에 의한 영향계수 Dij를 계산하기 위하여 꼭짓점을 공유하는 패널들이 이루는 각을 3차원 각도인 입체각을 이용하여 계산하였다.
입체각은 각을 3차원으로 확장한 것이며, 2차원 각이 평면에서의 퍼짐 정도를 나타낸다면, 입체각은 공간에서의 퍼짐 정도를 나타내는 척도가 된다. 입체각의 기호는 일반적으로 Ω로 나타내며, 라디안과 같이 무차원량이다. 입체각은 다음과 같이 표현 할 수 있다. 삼각형 패널의 꼭짓점을 Fig. 6처럼 벡터 R→1, R→2 그리고 R→3으로 정의하고, 삼각형 패널에 의한 O점에서의 입체각, Ω는 S/r2로 나타낼 수 있다. 여기서, S는 반지름이 r인 구위의 삼각형 ABC의 면적이고, r=1은 Ω=S임을 의미 한다.
삼각형 패널에 의한 입체각을 구하는 방법은 van Oosterom, A. and Strackee, J.(1983)에 의해 다음과 같이 정리 되었다.
여기서,
N=R→1∙(R→2×R→3)
D=R1R2R3+(R→1∙R→2)R3+(R→1∙R→3)R2+(R→2∙R→3)R1
자기유도 영향계수를 입체각으로 표현하기 위해 제어점과 특이점이 위치하는 꼭짓점을 포함하는 패널들이 이루는 입체각을 계산하게 된다. 제어점과 특이점을 포함하고 있는 꼭짓점에서 ε만큼 떨어진 점을 원점 O점으로 하여 각 패널별로 입체각을 구하고, 각 입체각을 더하여 해당 점에서의 자기유도 영향계수를 구한다.
Ωp=Ω1+Ω2+Ω3+Ω4+Ω5+Ω6=∑i6Ωi | (18)
|
여기서, i는 제어점과 특이점을 포함하고 있는 패널의 개수이다.
Fig. 7
Calculation of self-induction by solid angle
3.4. 압력 Kutta 조건
물체 주변의 순환(circulation)을 특정하기 위해 날개 뒷날 상, 하면에서의 압력이 같아지도록 압력 kutta 조건을 적용하였다. 비선형 압력 Kutta 조건은 Lee(1987)가 제안한 방법을 적용하였고 다음과 같다.
(∆ϕ)wake=Φu-Φl=ϕt.e.u-ϕt.e.l+U→∞ ∙r→t.e.+γ | (19)
|
여기서, ϕt.e.u와 ϕt.e.l는 날개 뒷날 상, 하면의 포텐셜 값이고, U→INF는 유입속도, r→t.e.는 날개 뒷날의 두 제어점 사이의 벡터, γ는 후류의 포텐셜 점프에 영향을 주는 계수이다.
후류의 포텐셜 점프에 영향을 주는 계수, γ는 날개 뒷날에서의 압력계수의 차이에 의해 결정되는 비선형항으로 반복계산을 통해 구할 수 있고, 본 연구에서는 secant 방법을 이용하여 구하였다.
γm(k)=γm(k-1)-∆Cpm(k-1)γm(k-1)-γm(k-2)∆Cpm(k-1)-∆Cpm(k-2), m=1…Mr | (20)
|
여기서, 첨자 k는 반복계산 횟수이고, ∆Cp는 날개 뒷날에서의 압력계수 차이다.
식 (20)을 식 (19)에 적용하면 다음과 같다.
(∆ϕ)m(k)=(ϕmNc+1-ϕm1)(k)+(U→∞∙r→t.e.)m+γm(k-1), m=1…Mr | (21)
|
이산화 식 (11)에 kutta 조건 식 (21)을 적용하면 다음과 같다.
∑j=1NnodeDijϕj(k)+∑m=1MrWim(ϕNc-ϕ1)m(k)=∑j=1NnodeSij(∂ϕ∂n)j-∑m=1MrWm[(U→∞∙r→t.e.)+γm(k-1)], i=1,2,...,Nnode | (22)
|