GPU 기반의 TSM-RT를 이용한 풍력 블레이드의 동적 RCS 계산

동적 RCS는 움직이는 표적의 시간에 따른 RCS의 변화를 의미한다^[12]. 움직이는 표적의 RCS는 식 (1)과 같이 표현할 수 있다.

여기서 t_p는 p번째 시간 샘플에서의 시간을, f는 입사파의 주파수를, 전체 장면을 구성하는 M은 산란점의 개수를 A_m(t_p)는 p번째 시간 샘플에서 m번째 산란점의 복소 진폭을, R_m(t_p)는 p번째 시간 샘플에서 m번째 산란점과 안테나 사이의 거리를, c=3×10⁸ m/s는 진공에서의 빛의 속력을 의미한다.

문제를 단순화하기 위하여 ‘레이다에서 방사된 전자파가 표적을 맞고 되돌아오는데 걸리는 시간은 같은 시간 동안의 표적의 움직임을 무시할 수 있을 정도로 충분히 짧다’고 가정한다. 따라서 각 시간 샘플 사이의 간격 Δt=t_p+1-t_p는 식 (2)와 같이 나타낼 수 있다.

여기서 R은 표적과 안테나 사이의 거리이고, PRF(Pulse Repetition Frequency)는 펄스반복주파수이다.

그림 1은 전자파 수치해석 방법을 이용한 동적 RCS 계산의 개념을 보여준다.

그림 1. | Fig. 1. 전자파 수치해석 방법을 이용한 동적 RCS 계산 방법의 개념 | Concept of dynamic RCS calculation using numerical method.

Download Original Figure

그림 1에서는 각 시간 샘플에 대한 3D 메쉬를 생성하고, 이에 대한 RCS를 계산하는 방법으로 동적 RCS를 획득한다. 이러한 개념으로 동적 RCS를 계산하기 위해서는 모든 시간 샘플에 대한 3D 메쉬가 필요하다. 그러나 1 kHz의 PRF를 가정할 경우, 1초의 동적 RCS를 얻기 위해 필요한 3D 메쉬의 개수는 1,000개에 달하므로 동적 RCS 계산에 필요한 3D 메쉬를 사용자가 모두 작성하는 것은 어렵다. 또한 1초의 동적 RCS를 얻기 위하여 1,000회의 RCS 계산을 수행해야 하므로 RCS 계산에 필요한 시간이 극단적으로 길어짐을 알 수 있다.

블레이드의 움직임을 표현하기 위해서 블레이드를 강체(rigid body)로 가정한다. 블레이드의 운동을 표현하기 위해서는 블레이드의 초기 위치와 회전 운동에 의한 위치의 변화를 블레이드의 3차원 메쉬에 반영해야 한다. 이것은 각 블레이드의 메쉬에 회전변환과 이동변환을 적용함으로써 표현이 가능하다. 그림 2는 움직이는 3차원 메쉬를 생성하는 절차를 보여준다.

그림 2. | Fig. 2. 움직이는 3차원 메쉬를 생성하는 절차 | Procedure to create a moving 3D mesh.

Download Original Figure

메쉬의 회전 변환은 식 (3)과 같이 계산이 가능하다.

여기서 $\overrightarrow{V}\in {ℝ}^3$는 메쉬를 이루는 정점(vertex)의 위치벡터이고, $\overrightarrow{R}\in {ℝ}^{3\times 3}$은 3차원 공간에서의 오일러 회전 행렬이다^[13].

메쉬의 이동 변환은 식 (4)와 같이 계산이 가능하다.

여기서 $\overrightarrow{M}\in {ℝ}^3$은 메쉬의 이동을 위한 벡터이다.

메쉬의 회전 변환과 이동변환은 동일한 연산을 메쉬를 이루는 수 많은 정점에 대해 수행하는 것이다. GPU를 이용하여 병렬화 된 선형대수 계산은 이미 많은 라이브러리들이 개발되어 있기 때문에 이러한 라이브러리를 사용하는 것만으로도 간단하게 GPU 기반의 병렬화가 가능하다. 본 논문에서는 NVIDIA에서 제공하는 GPU 기반의 선형대수 계산 라이브러리인 CUBLAS를 사용하여 해당 연산들을 구현하였다^[14]. 앞에서 설명한 방법으로 움직이는 3차원 메쉬를 표현하기 위해서는 각 부분의 메쉬에 회전 변환과 이동 변환을 어떠한 순서로 적용하고, 전자파 수치해석을 수행할 것인지를 지시하는 스크립트(script)가 필요하며, 이러한 스크립트는 변수, 함수, 조건문, 반복문과 같은 최소한의 문법을 지원해야 한다. 본 논문에서는 수치해석 코드에 프로그래밍 언어인 Lua의 인터프리터를 내장(embedding)하여 이 문제를 해결하였다^[15]. 따라서 3차원 메쉬의 움직임은 표 1과 같이 Lua 문법으로 표현이 가능하다.

표 1과 같은 스크립트를 이용하여 그림 2와 같은 절차를 표현할 수 있으며, 이러한 반복문에 표 1의 11번 줄과 같이 RCS 해석 함수를 배치하여 해석 결과를 배열에 저장하는 방식으로 움직이는 표적의 표현과 이에 대한 동적 RCS의 계산을 모두 수행할 수 있다.

동적 RCS를 계산하기 위해서는 레이다의 PRF와 관측시간을 곱한 횟수만큼의 RCS 계산이 필요하다. 물리광학법은 높은 주파수에서 RCS를 해석하는 기법들 중에서도 계산속도가 매우 빠른 방법이다. 물리광학법에 의하여 평면파가 입사하는 표면의 전류는 식 (5)와 같이 나타낼 수 있다.

여기서 ${\overrightarrow{H}}_i\in {ℂ}^3$는 입사파의 자계이고, $\overrightarrow{n}\in {ℝ}^3$은 평면파가 입사하는 표면의 법선 벡터이다.

식 (5)에서 $\overrightarrow{J}\in {ℂ}^3$는 입사파가 직접 도달하는 전면(lit region)과 도달하지 못하는 은면(shadow region)으로 나누어 계산하게 된다. 따라서 물리광학법으로 전류를 계산하기 전 은면을 미리 제거하는 과정이 필요하다. 이러한 은면 제거 과정은 물리광학법에 의한 RCS 계산 과정에서 가장 많은 계산시간을 요구하는 과정이다. 따라서 물리광학법의 계산 속도를 향상시키기 위해서는 반드시 광선추적법의 계산 속도를 향상시켜야만 한다. 입사파가 정의되었을 때, 은면을 제거하는 방법은 내적에 의한 방법, Z-buffer 방법, 광선추적법 등 많은 방법이 알려져 있으며, 최근에는 광선추적법 기반의 방법이 많이 사용된다^[16]. 물리광학법의 구현에서 광선추적법은 메쉬를 구성하는 각각의 삼각형 패치에 입사하는 광선들이 다른 삼각형 패치에 가로막히지 않고 도달할 수 있는 경우, 해당 삼각형 패치를 전면으로 그렇지 않은 경우는 은면으로 처리하는 방식이다. 이것을 그냥 계산할 경우, 하나의 삼각형 패치가 은면인지를 판단하기 위하여 메쉬를 구성하는 모든 삼각형 패치에 대한 충돌 검사를 수행해야 하고, 이는 매우 많은 계산량을 요구한다. 이러한 문제점을 개선하기 위하여 메쉬가 존재하는 공간을 분할하여, 충돌 계산시 공간과의 충돌을 우선 계산하고, 해당 공간에 위치한 삼각형 패치들만을 충돌계산에 참가시키는 방식의 가속 방법을 사용한다. 이러한 가속법의 대표적인 예시는 octree, kd-tree, bvh-tree 등이 알려져 있다^[17]～^[19]. 또한 이러한 광선추적법을 GPU 기반으로 병렬화시킴으로써 계산 속도를 향상시키는 연구들 또한 활발하게 이루어졌다^[18],^[19]. 본 논문에서는 NVIDIA의 GPU 컴퓨팅 라이브러리인 CUDA를 사용하여 GPU 기반의 octree를 구현하였다.

GPU를 이용하여 octree 기반의 광선추적법을 구현할 때의 어려운 점은 포인터에 의한 메모리 연결 구조와 연결리스트로 구성된 octree를 GPU로 전달하는 것이 어렵다는 것이다. 이러한 문제를 해결하기 위하여 그림 3과 같이 CPU에서 octree를 생성한 뒤, 포인터 대신 배열의 인덱스를 이용하여 octree의 각 노드를 방문할 수 있도록 한 배열의 형태로 전환하였다. 이러한 방법을 사용하여 성능저하 없이 GPU에서 동작하는 octree 기반의 광선추적법을 구현할 수 있다.

그림 3. | Fig. 3. Octree를 배열 형태로 바꾸는 방법 | Method to convert an octree into an array form.

Download Original Figure

그림 3에서 octree의 각 노드는 DATA ARRAY 배열에 깊이 우선 탐색의 방문 순서대로 저장된다. 이때, 각 노드데이터가 시작되는 지점과 각 노드의 길이를 각각 INDEX ARRAY 배열과 DATA LENGTH ARRAY 배열에 저장한다. 이러한 정보들을 이용하여 DATA ARRAY로부터 원하는 노드의 정보를 추출할 수 있다.

물리광학법을 계산할 때, 사용되는 삼각형 패치 한 변의 길이는 λ/8가 권장된다. 이는 각 삼각형 패치에서 발생하는 전계를 적분하기 위하여 필요한 조건이다. 따라서 이 조건은 주파수에 의존한다. 그러나 광선추적법에서 필요한 삼각형 패치는 3차원 모델의 경계를 정확하게 표현하는 것이 주된 역할이다. 따라서 광선추적법에서 요구되는 삼각형 패치의 크기 조건은 표적의 형상의 복잡도와 표면의 곡률과 같은 요소에 영향을 받을 뿐 주파수와는 무관하다. 따라서 광선추적법에서 요구되는 조건의 삼각형 패치를 사용하여 octree를 생성하고, 생성된 octree를 사용하여 기존의 λ/8 조건을 만족하는 메쉬의 충돌 검사를 수행하는 방법으로 은면 처리에 필요한 연산량을 크게 줄일 수 있다. 이러한 방법을 Two scale model ray tracing (TSM-RT)이라 한다^[20]. 추가로 물리광학법에서 단일 반사만을 고려할 경우, 입사광선 $\overrightarrow{s}$와 삼각형 패치의 법선벡터 $\widehat{n}$이 $\overrightarrow{s}\cdot \widehat{n}\le 0$인 경우, octree를 구성하는 메쉬에서 제외해도 무방하다. 이러한 방법으로 TSM-RT 방법의 연산량을 감소시킬 수 있다. 본 논문에서는 이 두 가지 방법을 모두 적용하였다.

완전 도체를 가정할 경우, 앞에서 설명한 방법으로 삼각형 패치 표면의 전류를 구하고 나면 식 (6)과 같은 원전계에서의 방사 적분 공식을 사용하여 전계를 구할 수 있다.

여기서 ω=2πf, μ는 진공에서의 투자율, S_n은 n번째 삼각형 패치의 면적, $\overrightarrow{r\prime }$은 삼각형 패치의 무게중심의 위치벡터, $\widehat{r}$은 삼각형 패치에서 관측지점으로 향하는 방향벡터, N은 표적을 이루는 삼각형 패치의 총 개수이다.

식 (5)의 Σ항을 각 패치별로 GPU에서 계산하도록 구현함으로써 GPU 기반의 가속이 가능하다. 본 논문에서는 CUDA를 사용하여 구현하였다.

본 논문에서 제안한 방법의 계산 결과를 검증하기 위하여 제안한 방법으로 계산한 결과를 상용 소프트웨어인 FEKO의 물리광학법으로 계산한 결과와 비교하였다. 본 논문의 모든 시뮬레이션은 Intel(R) Xeon(R) CPU E5-26700 @ 2.6 GHz Dual core CPU와, 128 GB 메모리, NVIDIA GeForce GTX 650 Ti 그래픽 카드를 장착한 환경에서 수행되었다. 제안된 방법과 FEKO를 이용한 시뮬레이션은 모두 동일한 환경에서 수행되었다. 그림 4는 구에 대하여 주파수를 변화시켜 가며 RCS를 계산한 결과와 계산에 걸린 시간을 보여준다. 그림 4에서 구에 대한 RCS 계산 결과는 FEKO와 제안한 방법으로 개발한 코드 (DARSIM)가 완벽하게 일치한다. 또한 단순히 비교하였을 때, DARSIM이 FEKO에 비하여 256배 이상 빠른 계산 속도를 보여준다. 그러나 FEKO는 주파수 변경 시 은면 계산을 다시 수행하고, DARSIM은 이전에 계산한 은면 처리 정보를 활용하는 등 동작 방식의 차이가 존재하므로 이러한 비교는 계산 시간에 대한 공정한 비교로 보기 어렵다.

그림 4. | Fig. 4. 구에 대한 RCS 해석 결과 | RCS analysis results for the sphere(frequency: 1～50 GHz, 256 points, radius: 0.1 m, patch number: 654,610).

Download Original Figure

그림 5는 임의의 곡면을 포함한 반경 1.95 m의 날개를 가진 풍력 블레이드에 대하여 블레이드의 회전 방향으로 각도를 변화시켜 가며 RCS를 계산한 결과와 계산에 걸린 시간을 보여준다. 개발된 코드에 전통적인 octree를 제외한 가속 기법을 적용하지 않았을 경우에는 FEKO보다 더긴 계산시간을 필요로 함을 확인할 수 있다. 여기에 GPU 가속을 적용할 경우, CPU로 계산했을 경우에 비하여 21.7 %의 계산시간만을 필요로 하며, 여기에 TSM-RT까지 적용할 경우 CPU로 계산했을 경우에 비하여 1.5 %의 계산시간만을 필요로 함을 확인할 수 있다.

그림 5. | Fig. 5. 날개가 3개인 블레이드에 대한 RCS 해석 결과 | Results of RCS analysis for blades with three wings(frequency: 10 GHz, theta: 0～360 deg, phi: 0 deg, 361 points, radius: 1.95 m, patch number: 465,360).

Download Original Figure

본 절에서는 개발한 코드를 사용하여 서로 다른 위치에서 두 개의 대형 블레이드가 함께 회전하는 장면에 대한 동적 RCS를 계산하였다. 표 2는 시뮬레이션에 사용한 파라미터를 보여주며, 표 3은 시뮬레이션에 사용된 표적의 속성을 보여준다. 그림 6은 계산된 동적 RCS의 시간 거리 그래프를 보여주며, 그림 7은 계산된 동적 RCS의 스펙트로그램을 보여준다.

그림 6. | Fig. 6. 시뮬레이션 결과(시간-거리 그래프) | Simulation result(time-range graph).

Download Original Figure

그림 7. | Fig. 7. 시뮬레이션 결과(spectrogram) | Simulation result(spectrogram).

Download Original Figure

표 3. | Table 3. 시뮬레이션 표적들의 정보 | Information of the simulation targets.

Download Original Figure

시뮬레이션에 소요된 시간은 1,304.3초이다. 표 2에서 입사파는 x축 방향에서 입사하므로 두 표적의 거리는 각각 5 m와 10 m에 위치해야 한다. 블레이드는 z축을 축으로 회전하며 한번 회전하는 동안 다가오고 멀어지는 블레이드에 의하여 총 6회의 정반사가 일어난다. 그림 6의 시간 거리 프로파일에서 표적의 거리가 각각 이론값과 일치하는 5 m와 10 m에 위치함을 확인할 수 있으며, 각표적에서 정반사에 의해 6회의 강한 RCS가 관측되는 간격이 각 표적의 1회전 주기와 일치함을 확인할 수 있다. 그림 7의 스펙트로그램에서는 두 표적의 도플러 성분이 겹쳐진 영상이 나타남을 확인할 수 있다. 스펙트로그램에서는 동일한 3개의 패턴이 하나의 주기를 이루며, 이 주기 역시 그림 7의 결과와 이론값이 일치함을 확인할 수 있다. 또한 두 표적의 블레이드 끝부분에 의해 발생하는 이론적인 도플러 주파수는 식 (7)과 같이 구할 수 있다.

여기서 λ는 레이다의 송신파장이며, 24 GHz에서는 0.0125 m이다. v_tip은 블레이드 끝부분의 속력이다. 식 (7)에 표 3에서 설정된 속성을 대입하여 각각의 표적에 대한 이론적인 최대 도플러 주파수를 구하면 16.7 kHz, 8.37 kHz와 같이 얻을 수 있으며, 이는 그림 7에서 관찰되는 최대 도플러 주파수와 일치한다. 이러한 결과는 제안한 방법으로 구현한 해석 소프트웨어가 둘 이상의 서로 다른 회전 표적의 동적 RCS를 올바르게 계산할 수 있음을 보여준다.