로터시스템 해석을 위해 국부 좌표계 기반의 변분 회전을 도입하여 전산해석을 단순화시킨 혼합 변분법 기하학적 정밀 보 구조 해석과 헬리콥터 로터 동적 해석을 위한 Peters-He가 제시한 동적 유입류 및 ONERA-EDLIN 동적 실속 모델을 결합한 로터 연계 해석 프로그램을 개발하여 적용하였다.

로터 시스템의 보 해석을 구성하기 위하여 Fig. 1과 같은 좌표계를 구성하여 좌표계 변환 행렬을 정의할 수 있다. Fig. 1에서 i는 고정 허브좌표계, a는 회전 허브좌표계, b는 보의 변형 전 초기 형상이 정의된 좌표계, 그리고 B는 변형 후 좌표계이다. 혼합변분법으로 정의된 보 지배방정식에서 경계조건 구성을 통해 회전 조인트 등의 다물체 요소 확장이 용이하다. 이 논문에서는 무힌지형 로터 블레이드를 모사하기 위하여 회전 조인트 및 베어링의 감쇠 및 스프링 계수를 적용하였다.

Fig. 1 
Present beam frame definition

Peters의 비정상 양력선 이론은 B-좌표계를 기준으로 정의되며, 식 (3)과 같이 좌표계 변환을 통해 자유류 속도와 유입류 속도가 정의된다. 유입류는 Peters-He의 3차원 동적 유입류를 적용하였다.

회전 변환 행렬(3×3)은 Euler 각 Φ 및 지수맵 회전 정식화를 이용하여 식 (4)와 같이 정의되며, 모든 물리량은 회전 허브 전역(material) 좌표계에서 정의된다.

Peters 비정상 양력선 이론은 정적 공기력을 풍동시험이나 CFD 계산을 통해 구한 정적 공기력 계수 C81표를 참조하여 계산하며, 정적 항력 계수를 이용한 표면 점성 항력 보정을 수행한다. 비정상 공기력은 식 (5)와 같이 계산한다.

ONERA-EDLIN 동적 실속은 식 (6)와 같이 동적 실속 추가 계수를 정적 공력 계수에 더하여 계산한다. 식 (7)에서 양력, 항력 및 모멘트 방향으로 각각 상태변수 형태의 미분방정식을 통해 동적 실속 추가 계수를 획득한다. 식 (6)에서 f_rev는 후퇴면 역흐름에 대한 변수이며, 익형 상대속도에 따라 0 또는 1의 값을 가진다.

코드 구성은 Eigen3 C++ 라이브러리를 이용하여 밀집 행렬 연산을 수행하고 Intel MKL PARDISO를 이용하여 희박행렬 계산을 수행하였다.

UH-60A의 해석에는 플랩, 래그, 비틀림 방향의 조인트, 래그 감쇠기 및 스워시 플레이트 등 다물체 구조를 모두 구성하였으며, 다물체 동역학 해석 프로그램 DYMORE 2.0을 활용하였다. DYMORE 2.0에 포함된 기능인 Peters-He 동적 유입류 및 양력선 이론을 사용하여 공기력을 계산하였다. 모델링 된 로터는 Fig. 2에 나타내었다. 이 연구의 해석 대상인 HART-II 로터와 UH-60A 로터의 제원을 Table 1에 정리하였다.

Fig. 2 
DYMORE2.0 UH-60A modeling

임의의 로터에 대하여 회전좌표계에서의 개별 블레이드 조종값을 식 (8)과 같이 선정하면 다중블레이드좌표계 회전변환 관계에 의해 비회전좌표계에서 N/rev의 출력만을 얻을 수 있음을 알 수 있다⁽⁸⁾. 괄호안의 지수는 블레이드 번호, 괄호 없는 지수는 N/rev, (N±1)/rev 주파수를 나타내며 각각 cosine 및 sine의 기저 함수로 구분할 수 있다. 식 (8)에서 블레이드 간 위상차가 없는 콜렉티브 조종 값은 N/rev 주파수를 그대로 전달하며, ±2π(i-1)/N, (i=1,2,3,…N)의 위상을 가지는 조종값은 (N±1)/rev의 주파수를 N/rev으로 각각 전달한다. 이러한 제어 입력에 대하여 고정좌표계에서의 허브 하중 변화량 또한 cosine 및 sine의 기저 함수로 나타내면 식 (9)와 같이 나타낼 수 있다.

전달함수 식별을 위하여 식 (8)의 각 기저 함수의 크기 성분을 로터의 동적해석 수행시 충분한 회전이 이루어진 정상상태에서 τ_threshold 시간 이후에 계단 입력으로 부여한다.

허브 하중 변화량의 N/rev 성분의 하중을 추출하기 위하여 N/rev 주파수에 해당하는 신호의 극값들을 추출하여 계단 응답에 대한 전달함수 추정을 수행하였다. 식 (11)과 같이 허브 수직전단 하중과 롤링, 피칭 모멘트에 해당하는 6개 성분만을 사용하여 전달함수는 6개 입력 벡터 성분 θcN,θsN,θcN+1,θsN+1,θcN-1,θsN-1T과 6개 출력 하중 δFzcN,δFzsN,δMxcN,δMxsN,δMycN,δMysNT으로 구성하였다. 각 입력성분에 해당하는 SISO 전달함수는 zero를 포함하는 6차/6차 형태로 모델링이 가능하다. 전달함수를 대각행렬로 정의하고 이를 상태공간 방정식으로 변환하여 LTI 시스템 행렬 A, B, C, D를 획득할 수 있다.

식별된 선형 시스템 행렬 A, B, C, D를 이용하여 LQG(Linear quadratic Gaussian) 제어 구조를 적용하면 계측되는 N/rev 로터 진동 하중에 대해 N/rev 섭동 하중을 발생시켜 로터 진동을 제어할 수 있다.

Fig. 3은 로터 진동을 감소시키기 위한 블록선도를 나타낸다. MATLAB/Simulink 환경에서 시뮬레이션을 수행하여 정적 게인을 획득하고 이를 해석 프로그램에 적용하여 제어 해석을 수행하였다. 상태변수의 관측기는 정상상태 칼만 필터를 적용하였다. 관측기 이득은 트림 해석결과의 FzcN,δFzsN,MxcN,MxsN,MycN,MysNT의 N/rev 하중성분의 크기를 대각행렬 R_k의 구성성분으로 하고, 대각행렬 Q_k의 크기를 반복적으로 관측기의 성능을 만족하도록 구할 수 있다. 제어 이득을 획득하기 위하여 식 (12)와 같은 출력 궤환 방식의 목적함수를 이용하여 Q와 R을 선정하였다. 식 (12)에서 HL_i는 i-방향의 허브 진동 하중이며 |u_max|는 IBC 제어 입력 최대 크기이다. α_i, β_i및 ρ는 이 연구에서는 최적화를 수행하지 않았으며 휴리스틱 방법으로 결정하였다.

Fig. 3 
LQG control diagram for the present IBC system

Fig. 4는⁽⁹⁾ HART-II의 5.3° 하강비행 조건(baseline)에서 로터 블레이드의 끝단 변형 크기를 나타내며, Fig. 5는⁽⁹⁾ 블레이드의 뿌리로부터 17 % 지점에서 플랩 방향 굽힘 모멘트를 나타낸다. 개발된 해석 프로그램이 시험 결과와 유사하며, 동일한 수준의 보 이론 및 공기역학 이론을 사용한 기존 로터 통합 해석 프로그램과 비교하여 유사한 수준에서 예측하고 있음을 발견하여 개발된 해석 프로그램의 유효성을 확인하였다⁽⁹⁾.

Fig. 4 
Tip flap-wise displacement of the baseline HART II

Fig. 5 
Flap-wise bending moment at 0.17R of the baseline HART II

Fig. 6은⁽⁹⁾ HART-II 로터의 회전속도에 따른 모드 해석 결과를 나타낸다. 상용 해석 프로그램과 유사한 예측 성능을 보이고 있으며, 동일한 요소 개수에서 요소간 물성 보간에 선형 보간을 사용하는 CAMRAD-II, DYMORE 4.0 보다 더 실험 결과에 가까운 예측 결과를 보였다.

Fig. 6 
HART II fan plot

UH-60A 로터는 비행 시험 결과의 C8513 조건에서 해석을 수행하였다. 이는 전진비 0.15, 로터 디스크 하중 계수 0.08 및 0.76° 하강비행 조건으로 HART-II의 해석 조건과 유사하다. Fig. 7 (a) ~ (d)에서 이 연구에서 구성된 UH-60A 로터의 피치링크 강성을 조절함에 따라 참고문헌들의 결과를 재생산할 수 있음을 보이며, Fig. 8 (a), (b)는 블레이드 30 % 지점과 50 % 지점에서 트림 수행후 플랩 방향 굽힘 모멘트를 나타낸다^(10~13). 따라서 현재 사용된 DYMORE 2.0 해석파일이 올바르게 구성되었음을 관측하였다.

Fig. 7 
UH-60A fan plot

Fig. 8 
Flap-wise bending moment of the baseline UH-60A

Fig. 9 (a)와 (b)에 HART-II 로터와 UH-60A 로터에 대해 도출된 전달함수 결과 일부를 도시하였다. N/rev 주파수 성분의 과도응답의 특성이 적절히 모사되었음을 확인할 수 있고, 1.5/rev 이하의 주파수 동특성 범위를 가진다.

Fig. 9 
Transfer function identification results

Fig. 10은 HART-II 로터의 제어 해석 결과를 나타낸다. 로터의 트림 해석과 진동하중 제어기가 코드 내에서 동시에 작동하도록 하였으며, LQG 제어기가 작동하기 이전에 하중의 변화 구간은 로터 트림을 위한 콜렉티브 및 사이클릭 피치의 민감도 행렬을 식별하는 구간이다. LQG 제어기 작동시 추력이 일부 감소하나 트림제어기의 작동으로 트림상태로 복귀함을 보이며, 트림 제어기가 작동할 경우에도 평균 50 %의 진동하중 감소 성능을 보였다.

Fig. 10 
HART-II vibratory load control results

Fig. 11은 UH-60A의 진동하중 제어 예측 응답을 도시하고 있다. DYMORE2.0을 활용하여 MATLAB/Simulink 선형 시스템에서 도출된 제어 입력을 개루프 방식으로 인가하는 방식으로 해석을 수행하였다. 롤링, 피칭 모멘트 방향의 모델링 오차가 다소 존재함을 관측하였다.

Fig. 11 
UH-60A vibratory load control results

선형 모델링된 시스템에서 예측한 성능과 비선형 시뮬레이션에 적용된 제어기의 응답 예측이 유사한 수준임을 확인하였으며, 선형시스템의 위상 및 이득 여유가 적절히 설계될 수 있음을 의미한다. MATLAB/Simulink의 루프 절단 방법을 통해 개루프 시스템 전달함수를 얻어 다중 입출력 시스템에 대한 상대 안정도를 디스크(disk) 여유 계산을 이용하여 해석하였다. 계산된 각 로터의 최소 여유값은 Table 2에 요약하였다.