n 자유도 점성감쇠 시스템의 운동방정식은 아래와 같다.

여기서, M은 질량행렬, C는 감쇠행렬, K는 강성행렬이며, x(t)와 f(t)는 시간에 따른 변위 벡터와 힘 벡터이다. 위첨자 ‘C’는 감쇠시스템을 의미한다. 조화가진, f(t)=F(ω^jwt),으로 가정하면 식 (1)은 아래와 같이 정리된다.

여기서, x(t)^C=X(ω)^Ce^jwt이며 j는 허수를 나타낸다. 식 (2)에서 K-ω²M은 아래와 같이 표현될 수 있다.

여기서, α(ω)^N은 실수 혹은 정상상태 주파수 응답함수을 나타내며, 위첨자 ‘N’은 실수 혹은 정상상태를 의미한다. 식 (3)을 식 (2)에 대입하고, 양변에 [α(ω)^N]를 곱하여 식 (4)를 얻을 수 있으며, 그 후 식 (5) 변환행렬[9]를 식 (4)에 대입하여 정리하면 식 (6)을 얻을 수 있다.

식 (6)에서 변위 벡터, X(ω)^C,는 힘 벡터와 복소 주파수 응답함수의 곱으로 표시될 수 있다.

식 (7)을 이용하면, 식 (6)은 아래와 같이 정리된다.

식 (8)을 실수부와 허수부로 정리하면

정상상태 주파수 응답함수의 경우 허수부는 0이 되어야 하므로,

식 (10)을 감쇠행렬로 정리하면,

식 (11)은 감쇠행렬을 구하기 위한 기본식이 되며, 정상상태 주파수 응답함수와 복소 주파수 응답함수의 관계는 다음과 같다.

측정을 통해 m개의 복소 주파수 응답데이터를 획득할 경우 식 (11)은 아래와 같이 정리될 수 있다.

여기서, 위첨자 ‘+’는 유사역행렬(Pseudo inverse matrix)을 의미한다. 정리하면, 측정을 통해 m개의 복소 주파수 응답데이터를 n개의 자유도에서 획득할 경우 식 (13)을 이용하여 [n×n] 크기의 감쇠행렬을 추정할 있다.

제안된 기법을 검증하기 위하여 m₁,m₂,m₃는 10kg, 14kg, 12kg이며, k₁,k₂,k₃는 2000N/m, 3000N/m, 2500N/m이고, c₁,c₂,c₃는 2.1Ns/m, 3.2Ns/m, 2.5Ns/m인 3자유도 비비례 점성감쇠 시스템을 Figure 1과 같이 구성하였다. 이 시스템을 구성하는 M, C, K는 Table 1과 같다. 이 중 M과 K는 알고 있는 상태로 가정하였으며, C는 모르는 상태로 제안된 기법을 통해 추정하는 수치검증 시나리오를 고려하였다. 3자유도 수치모델로부터 9개의 복소 주파수 응답함수를 획득 하였으며, 이중 0 ~ 50Hz 범위에서 획득된 복소 주파수 응답함수 데이터만을 감쇠행렬 추정에 이용하였다.

Figure 1: 
3-DOF nonproportional viscous damping system

식 (13)을 통해 추정된 시스템의 복소 주파수 응답함수와 수치모델의 복소 주파수 응답함수를 비교하여 정확성을 검토하였다. Figure 2와 Figure 3은 추정된 모델과 수치모델의 복소 주파수 응답함수, FRF_1,1과 FRF_2,3 이다. FRF_i,k에서 i는 가진점 위치, k는 응답점 위치를 나타낸다. 검정색 실선은 수치모델의 복소 주파수 응답함수이며, 빨간색 점선은 추정된 시스템의 복소 주파수 응답함수로 두 선이 매우 잘 일치하는 것을 확인 할 수 있다. 추정된 감쇠행렬과 오차율을 Table 2에 정리하였다. Table 2로부터 추정된 감쇠행렬이 수치모델과 정확히 일치하는 것을 확인 할 수 있다.

Figure 2: 
Overlay plot of simulated experimental FRF_1,1 and identified FRF_1,1

Figure 3: 
Overlay plot of simulated experimental FRF_2,3 and identified FRF_2,3

제안된 기제안된 감쇠행렬 추정기법을 검증하기 위하여 원형 외팔보 모델을 Figure 4와 같이 생성하였다. 수치해석 모델에 사용된 물성치와 수치는 Table 3과 같다. 총 20개의 요소와 21개의 절점으로 구성되었으며, 각 절점 당 2개의 자유도(회전과 병진 자유도)를 가지고 있다. 따라서 한쪽 고정조건을 고려하면 외팔보의 질량행렬과 강성행렬은 38×38의 크기를 가진다. 이러한 경우 식 (13)을 통해 감쇠행렬을 추정하기 위해서는 총 1444개의 복소 주파수 응답함수를 추정해야 하는 어려움이 있다. 또한 실제 실험모델의 자유도는 센서의 측정자유도와 사용된 개수의 곱으로 제한되기 때문에 자유도 축소가 고려되어야 한다. 이러한 자유도의 제한을 고려하기 위해 자유도 축소기법[10]을 사용하여 3, 5, 7, 9, 11, 13, 15, 17, 19, 21번 절점에서 각각 3, 7, 11, 15, 19, 23, 27, 31, 35, 39번의 병진자유도만을 가지는 10자유도 모델로 수치모델을 축소시켰으며, 축소된 모델과 축소 전 모델의 1 ~ 10차 고유주파수들의 차이가 0.05 %보다 작게 되도록 하였다. 축소된 질량행렬과 강성행렬을 각각 M_R과 K_R로 정의하였으며, 식 (14)를 이용하여 비비례 점성감쇠 행렬을 구성하였다.

Figure 4: 
Finite element model of a cantilever beam

여기서, C_i,j와 K_Ri,j는 감쇠행렬과 강성행렬의 i번째 행과 j번째 열의 성분, α_i,j는 0 ~ 1사이의 불규칙 값이다. 생성된 비비례 감쇠행렬을 Table 4에 나타내었다.

M_R, K_R, C로부터 100개의 복소 주파수 응답함수를 획득하였다. 이중 0 ~ 5000Hz 범위에서 획득된 복소 주파수 응답함수 데이터(Case 1)와 0 ~ 1000Hz범위에서 획득된 복소 주파수 응답함수 데이터(Case 2)를 각각 감쇠행렬 추정에 이용하였다. Case 1과 Case 2에서 추정된 모델과 수치모델의 복소 주파수 응답함수를 Figure 5와 Figure 6에서 비교하였다. Case 1과 Case 2 모두에서 추정된 모델의 복소 주파수 응답함수와 수치모델의 복소 주파수 응답함수가 매우 잘 일치하는 것을 확인 할 수 있다. 이 결과는 직접법이 가지고 있는 장점을 잘 보여준다. [n×n] 크기의 감쇠행렬을 추정하기 위해서 직접법은 [n×n×m] 크기의 복소 주파수 응답함수 행렬이 필요하고, 간접법의 경우 n개의 모드 특성치가 필요하다. 다시 정리하면, [n×n] 크기의 감쇠행렬을 추정하기 위해서 직접법의 경우 [n×n] 점에서 복소 주파수 응답함수를 계측하면 되기 때문에 간접법과는 다르게 자유도에 제한이 없다. 간접법의 경우 일반적으로 [n×n] 크기의 감쇠행렬을 추정하기 위해서는 n개의 고유주파수 혹은 고유모드 등을 추정하여야 하는 단점이 있다.

Figure 5: 
Overlay plot of simulated experimental FRF_1,1 and identified FRF_1,1,, Case 1

Figure 6: 
Overlay plot of simulated experimental FRF_1,1 and identified FRF_1,1,, Case 2

진동실험을 통해 실제 계측된 데이터에는 잡음이 필연적으로 포함되기 때문에 잡음에 대한 강건성(Robustness) 검토가 필요하다. 잡음에 대한 강건성 검토를 위해 Case 1의 복소 주파수 응답함수 데이터에 2%의 잡음을 인위적으로 추가하여(Case 3) 감쇠행렬 추정에 적용하였다. Figure 7은 Case 3에서 추정된 모델의 복소 주파수 응답함수와 수치모델의 복소 주파수 응답함수를 비교한 것으로 매우 잘 맞는 것으로 보인다. 그러나 Figure 7에서 보이는 것과는 달리 추정된 감쇠행렬(Table 5)은 수치모델의 감쇠행렬과 많은 차이를 보였는데, 이러한 오차는 복소 주파수 응답함수에 포함된 잡음에서 기인한 것으로 생각된다.

Figure 7: 
Overlay plot of simulated experimental FRF_1,1 and identified FRF_1,1, Case 3