Modeling and Bifurcation Analysis of the 2D Airfoil with Torsional Nonlinearity

비틀림 비선형성을 갖는 2차원 익형의 모델링 및 Bifurcation 해석

Abstract

Recent developments for high altitude, long endurance conventional UAVs(HALE UAVs) have revealed new issues regarding aircraft structure design and analysis. First of all, due to intensive mission requirements, the structures of HALE UAVs have lightweight and very flexible main wing with high aspect ratio, and slender fuselage. For this kind of structures, aeroelastic characteristics are different from conventional aircrafts. Hence, currently developed analysis methods are not suitable to fully understand strucutral dynamics of the very flexible aircraft, and to guarantee structural reliability. Therefore, various structural studies considering nonlinear behaviors which are generally ignored for the conventional aircraft strucutral analyis have been attracting researchers interests. Nonlinear flutter of the very flexible wing is one of the subject to be studied in combination with strong coupling between aeroelastic characteristics and flight dynamics. Herein, as preliminary study, modeling and nonlinear system analysis of the 2D airfoild with torsional nonlinearity have been discussed.

Nomenclature

h Plunge(플런지) 자유도 α Pitch(피치) 자유도 m 단위 길이당 질량 Iα 비틀림 관성 모멘트 p 공기 밀도 kh 선형 굽힘 스프링 상수 ka 선형 비틀림 스프링 상수 Gh 비선형 굽힘 스프링 상수 Ga 비선형 비틀림 스프링 상수 U 유속

 

1. 서 론

최근 항공 선진국들 주도로 개발되고 있는 고고도 장기 체공 무인기의 강화된 임무 수행 조건을 만족시키기 위해, 보다 높은 종횡비의 날개와 가늘고 작은 동체형상을 갖는 항공기 구조 설계의 필요성이 증가하고 있다.

2003년 미국에서 행해졌던 Helios의 항공시험 중 파손 사건은(1), 유연항공기 구조설계의 신뢰성 확보를 위한 새로운 구조해석 방법의 연구를 촉진시키는 계기가 되었다. 이렇듯, 높은 종횡비를 가진 경량, 유연구조 비행체의 공탄성 특성은 기존의 선형해석 만으로는 해석이 불가능하게 만드는 요인들(large deformation, low mode frequency 등)을 포함하고 있으며, 이를 극복하기 위한 비선형 해석 방법들이 연구되어지고 있다(2,3). 최근에는 유연항공기 주익의 Hopf bifurcation에 의한 limit cycle flutter에 대한 연구 또한 진행되고 있다(4,5).

이 논문에서는 limit cycle flutter 해석에 널리 사용되는 단순한 2차원 익형 모델(6~9)을 활용하여 구조적 비틀림 비선형성을 포함한 2차원 익형의 모델링 및 Hopf bifurcation에 의한 비선형 플러터 거동에 대한 비선형 시스템 해석을 다루었다.

 

2. 모델링

2.1 개념도 및 운동방정식

Fig. 1에 2차원 익형 모델의 개념도를 제시하였다. 이 논문에서 다루고자 하는 2차원 익형은 조정면을 포함하고 있지 않으며 plunge (h)와 pitch (α)의 2자유도를 갖는다. 여기서 L 과 M은 각각 날개에 걸리는 양력(lift)과 비틀림 모멘트(torsional moment)이며, unsteady 유동의 양력은,

Fig. 1Schematics of 2D airfoil

으로 나타낼 수 있으며, 비틀림 모멘트는,

와 같이 표현할 수 있다(10). 식 (1), (2)에 있는 C(k)는 Theordoson's function이며, 이 논문에서는 quasi-steady 유동을 가정하여 C(k)=1로 정의한다. 공탄성 거동의 비선형 모델링을 위해 궁극적으로는 plunge와 pitch 모두의 비선형 성을 가정하는 것이 맞지만, 기초연구를 위한 개념 모델로서는 pitch 비선형성만을 가정한 모델들 또한 여러 연구에서 많이 사용되고 있으며, 이 해석 모델에서도 모델링, 무차원화 및 비선형 거동의 단순화를 통한 해석의 용이함을 위하여 pitch의 3차항에 의한 hardening stiffness만을 가정하도록 한다. 이를 반영한 plunge와 pitch의 spring force는 각각,

으로 표현할 수 있다. 위에서 설명한 바와 같이 해석 모델의 비틀림 비선형성과 quasi-steady airflow 상태를 가정하였을 때, 2차원 익형 모델의 최종 운동 방정식은,

와 같이 나타낼 수 있다(11).

2.2 무차원화 및 상태방정식

수치해석을 위해, 운동방정식을 무차원화 시키고 그 후 상태방정식으로 표현하는 과정을 살펴보면 다음과 같다. 무차원화된 plunge, pitch, time은 각각 , α, 로 정의한다. 이를 이용하여 식 (4) 를 무차원화 하면,

와 같다. 우측 항의 과 은 각각 무차원화된 양력과 비틀림 모멘트이며,

와 같이 표현할 수 있다. 여기서 정의된 무차원 상수들은 Table 1에서 설명하였다. 이후의 선형 및 비선형 해석과 수치해석을 위해, 식 (5)와 (6)을 상태 방정식으로 표현하는 것이 용이하다. 우선, 다음과 같이 F1과 F2을 정의한다.

Table 1Nondimensional parameters

2차 미분항을 한 개 씩만 포함한 두 식으로 분리하기 위해 식 (7)의 각 함수에서 와 를 각각 소거 하면 다음과 같이,

와 를 각각 포함한 Fα 두 개의 식으로 분리할 수 있다. 식 전체의 전개는 이 논문상에는 생략 하도록 한다. 총 4개의 상태 변수를 {𝓎l,𝓎2,𝓎3,𝓎4}T = 와 같이 정의하고, 식(8)의 두 식을 이용하여 상태 방정식을 구성하면,

와 같이 표현할 수 있다. 여기서 행렬 [J] 는 상태변수의 1차항의 계수들로 이루어진 행렬이며, 안정평형점인 원점 ({𝓎l, 𝓎2, 𝓎3,𝓎4}T={0,0,0,0})에서의 시스템의 Jacobian이다. 벡터 {N}은 pitch non linearity 의 계수이다. Table 1에 제시된 무차원 상수들로 이루어진, [J] 와 {N}에 나타난 각 계수들을 간단한 알파벳으로 표현하면 다음의 구조를 나타내고 있다. [J] 는

와 같이 표현되며, {N}은

와 같이 표현할 수 있다. 문자로 표현되어진 각 계 수항들은 부록에 수록하였다.

 

3. Bifurcation 해석

3.1 Bifurcation Point

시스템의 평형점이 비선형 불안정 상태로 전환되는 해당 상수(제어상수)의 임계값을 bifurcation point라고 한다. Bifurcation point, 즉 이 모델에서의 플러터 점을 구하기 위해서 평형점의 Jacobian을 먼저 구하여야 한다. 위에서 정의한 상태공간에서의 평형점은 원점이기 때문에 이는 자연스럽게 식 (9)에서의 [J] 가 된다. 평형점의 안정성 판별은 평형점 부근의 선형화 영역에서 이루어지므로, bifurcation point를 구할 때 비선형 항들은 고려하지 않는다. 원점 y0={0,0,0,0}T 의 안정성 판별을 위한 특성 방정식은,

로 표현할 수 있다. 이는 Jacobian ( )의 eigenvalue를 계산하여 평형점의 안정성(stability)을 판별하는 것으로, complex conjugate pair인 두개 의 eigenvalue가 허수축을 좌측에서 우측으로 가로지를 때(실수부가 0이 될 때, 즉 λ1,2± jw) 평형점에서 Hopf bifurcation이 일어나게 된다. 이 논문에서는 제어상수로 V(nondeimensional freestreem velocity), 𝑥α(distance from elastic axis to c.g.), 그리고 μ(density ratio)의 세 가지를 고려하였다. 제어 상수를 제외한 시스템 상수들은 각각, rα=0.5, =0.5, Gα=0.5, α=−0.35를 사용하였다. Fig. 2는 각 제어상수들의 변화에 따른 특성방정식의 근궤적 선도를 보여주고 있다. 각각의 근궤적 선도에서 첫번째 complex conjugate eigenvalue가 허수축을 가로지를 때 Hopf bifurcation이 일어난다. 이처럼 단 순화된 모델에 대해서는 직접 Jacobian을 구하고, eigenvalue들을 추적함으로써 bifurcation point를 구 하는 것이 가능하다.

Fig. 2Root locus plot for each control parameter

Hopf bifurcation은 supercritical Hopf bifurcation과 subcritical Hopf bifurcation으로 분류할 수 있으며, 이 논문에서 처럼 3차항에 의한 hardening stiffness를 가정하였을 경우 supercritical Hopf bifurcation( Fig. 3)이 일어나게 된다.

Fig. 3Supercritical Hopf bifurcation

3.2 Bifurcation Boundary

위에서 살펴보았던 세 가지의 제어 상수 V, 𝑥α,μ에 대하여 parameter space상에 시스템의 안정 영영과 불안정 영역을 나타낼 수 있다면, 시스템 상수들을 운용 목적에 맞게 최적화 시킬 수 있는 길잡이가 될 수 있다. 이처럼 parameter space상에서 시스템의 안정/불안정 영역을 나누는, bifurcation point들로 이루어진 경계선을 bifurcation boundary라 한다. Bifurcation boundary는 시스템 모니터링 및 센서 등의 다양한 활용에 연구되고 있으며(12), 이 논문에서는 시스템의 파라미터 개념설계의 기법으로서 접근해 보도록 한다.

Bifurcation boundary는 이론식, 수치해석, 실험단계에서 모두 가능하다. 실제 존재하는 시스템에 대하여 제어가 가능한 상수들로 이루어진 parameter space상에 실험을 통해 구해진 bifurcation point들을 표시하여 선도를 그릴 수 있으며, 복잡한 모델식에 대해서는 각각의 상수에 대하여 numerical integration을 통해 직접 bifurcation point를 찾아 실험과 같은 방법으로 bifurcation boundary를 얻을 수 있다.

이 논문에서는 단순화 된 모델의 이론식으로부터 bifurcation boundary를 찾는 과정을 다루었다. 앞에서도 논의하였듯이, 식 (12)의 특성방정식으로부터 시스템에 Hopf bifurcation이 일어나는 조건은 한쌍의 complex conjugate eigenvalue가 복소수평면의 허수축상에 위치할 때, 즉 λ=± jw일 때 이다. 이를 특성방정식에 대입하면 다음과 같이 두 개의 식을 얻게 되며,

현재 시스템상에서 아직 정의되지 않은 값들은 V, 𝑥α,μ, 그리고 bifurcation이 일어나는 frequency w이다. 𝑥α−V parameter space상에서 bifurcation boundary를 구하기 위해 μ를 특정 값으로 정의한다. 여기서는 μ=7, μ=9, 그리고 μ=11의 세 가지 값을 고려하도록 한다. μ를 결정한 후, w를 0부터 충분히 작은 𝞓w 간격으로 증가시켜가며 각 w값에 대하여 식 (13)으로 부터 𝑥α, V 값을 계산한다. w는 운용범위 내의 parameter space 상에 boundary가 형성되는 값들을 그 범위로 지정한다. 주의해야 할 점은, 𝑥α, V 값이 계산되었을 때 ± jw이외의 eigenvalue들이 모두 복소수 평면의 허수축 좌측에 위치하고 있는지 확인해 보아야 한다는 것이다. 만일 허수축 우측에 다른 eigenvalue들이 존재하고 있다면 시스템은 이미 불안정 상태이며, 이때 주어진 w 와 구해진 해 𝑥α, V 는 모두 bifurcation과는 무관한 값이다.

위의 과정을 통해 얻은 bifurcation boundary는 Fig. 4에서 확인할 수 있으며, 이를 통해 limit cycle flutter가 일어나는 조건과 고려된 제어 상수들의 상관관계를 경계선도로 나타낼 수 있다. Fig. 4를 살펴보면, 무게중심이 elastic axis에 가까워질수록, 즉 𝑥α가 감소할수록 안정상태의 유속 범위가 커지고, 날개의 질량이 증가할수록, 즉 μ가 증가할수록 안정상태의 유속 범위가 커진 다는 것을 쉽게 알 수 있다. 이처럼 개념 설계 단계에서 bifurcation boundary를 통해 목표로 하는 운용 범위(유속 범위) 내에서 시스템 상수들(질량, 무게중심의 위치 등)의 1차적인 최적화를 수행 할 수 있다.

Fig. 4Bifurcation boundary on 𝑥α−V parameter space

 

4. 결 론

간단한 수학적 모델을 통해 알 수 있듯이, 비선형 항은 시스템의 플러터 점, 즉 bifurcation point의 위치에 영향을 끼치지 않으나, 시스템에 bifurcation이 일어난 이후의 비선형 거동의 경향에는 직접적인 영향을 끼친다. 이 논문에서처럼 3차항에 의한 hardening stiffness를 고려하였을 경우 supercritical Hopf bifurcation이 일어나게 된다.

2차원 모델을 통해 제어 상수들을 정의하고 parameter space 상에서 bifurcation boundary를 얻어냄으로써 상수들에 대한 1차적인 최적화가 가능하다. 물론, 이 연구는 모델 자체의 단순성으로 인해 실제 유연항공기의 구조 비선형성으로의 직접적인 적용은 적합하지 않은 한계를 갖는다. 그러나, 단순화된 모델을 바탕으로 추후 비선형성이 고려된 2차원 익형의 풍동 시험을 수행하고, Hopf bifurcation의 예측 인자(13) 등의 최신 비선형 시스템 해석 이론을 플러터 현상에 적용하는 등의 추가적인 이론 연구를 수행할 수 있는 기반으로서의 의미가 있다고 할 수 있다.