구조물의 모드 계산을 위해 고려되는 비감쇠 자유진동 운동방정식은 다음과 같다.

여기서 M_g과 K_g는 각각 질량 및 강성행렬을 나타내며, u¨g와 u_g는 각각 가속도 및 변위 벡터를 나타낸다.

식 (8)의 해인 변위 벡터 u_g 와 가속도 벡터 u¨g는 다음과 같은 형태로 가정할 수 있다.

여기서 x_g는 조화 운동 sin (ωt+θ)에 대한 크기 벡터(magnitude vector)를 의미한다.

식 (9)와 식 (10)을 식 (8)에 대입하면 다음과 같이 정리된다.

식 (12)은 고유치 문제(eigenvalue problem)와 동일한 식이므로 다음식과 같이 재표현 할 수 있다.

여기서, λ_i와 (ϕ_g)_i는 각각 i 번째 고유값(eigenvalue)과 고유벡터(eigenvector) 이고, N_g는 고려한 모델의 총 자유도 개수이다. 일반적으로 고유벡터(eigenvector) (ϕ_g)_i를 고유 모드(eigen mode)라 표현하기도 한다.

앞서 계산된 고유 벡터(eigenvector)들을 사용하여, 다음 식 (14)와 같이 변위 벡터 u_g를 정의할 수 있다.

u_g = Φ_gq_u with

여기서 Φ_g는 고유벡터(eigenvector) (ϕ_g)_i들을 포함하는 고유벡터 행렬(eigenvector matrix)이고, q_u는 모드 기여 계수 벡터를 나타내며, 일반좌표벡터(generalized coordinate vector)라 하기도 한다.

식 (14)로부터 주지해야 할 사실은 고유벡터(eigenvector) (ϕ_g)_i들은 변위 벡터 u_g를 표현하기 위한 기저 벡터(basis vector)로 사용된다는 점이다.

식 (14)에 무어-펜로즈 유사 역행렬(moore-penrose inverse)을 사용하여 일반좌표벡터 q_u는 다음과 같이 표현된다.

위첨자 *는 무어-펜로즈 역행렬을 의미하며, 일반적으로 무어-펜로즈 역행렬은 직사각형 행렬의 역행렬을 수행하기 위해 사용된다.

식 (14)에서 표현된 바와 유사하게 (즉, ϕ_g = Φ_gq_u), 변형률 벡터 ϵ_g 를 표현하기 위한 기저 벡터(Basis vector)가 존재하고, 그에 상응하는 기여 계수가 존재한다.

ϵ_g = Ψ_gq_s with

여기서, Ψ_g는 변형률 모드 행렬(strain mode matrix)이고, 그 안에 포함되어 있는 벡터 (ψ_g)_i 들은 변형률 모드(stain mode)라 한다. 이 벡터들은 변형률 벡터 ϵ_g를 표현하기 위한 기저 벡터로 사용된다. 그리고 q_s는 각 변형률 응답에 기여하는 변형률 모드 기여 계수(modal participation factor) 벡터이다.

변위 모드 기여계수 q_u와 마찬가지로, 식 (16)에 변형률 모드 기여 계수 q_s는 무어-펜로즈 유사 역행렬을 적용하여 다음 식 (17)과 같이 표현할 수 있다

또한, 변위 모드 기여계수 q_u와 변형률 모드 기여 계수 q_s는 동일하다고 가정하면, 다음 식 (18)과 같이 표현할 수 있다.

따라서, 식 (15), (16), (17), (18)로부터, 다음과 같이 변형률 벡터 ϵ_g와 변위 벡터 u_g에 대한 관계식을 도출할 수 있다.

최종적으로, 식 (19)로부터, 다음과 같이 변형률 벡터 ϵ_g와 변위 벡터 u_g에 대한 관계식 (20)을 얻을 수 있다.

여기서 행렬 B는 일종의 계수 행렬(coefficient matrix)이며, 변형률 벡터 ϵ_g와 변위 벡터 u_g 사이의 관계를 정의할 수 있는 변형률-변위 관계 행렬(strain-displacement relation matrix)라 정의할 수 있다.

시간 t+Δt에서 평형 상태의 구조물의 동적 방정식은 다음 식 (21)과 같이 표현된다.

여기서, C_g는 감쇠 행렬(Damping matrix)를 나타내고, F는 시간 하중 벡터를 나타낸다. 이 식의 해는 직접시간적분법(direct time-integration method)을 통해 도출할 수 있으며, 본 연구에서는 Newmark-β method (Newmark, 1959)을 적용하여 해를 도출하였다.

Newmark-β method로부터 도출한 해인 가속도 벡터 u¨g와 속도 벡터 u˙g는 다음과 같이 표현된다.

여기서 β와 γ는 원하는 정확도와 안정성에 따라 결정될 수 있는 매개변수이다. Newmark는 인위적인 감쇠를 피하기 위해 γ = 0.5 값을 제안했다. β의 값은 시간 간격 t와 t+Δt 동안 가속도가 변한다고 가정되는 방식에 따라 변하고, Trapezoidal rule (Newmark, 1959)에 의해 β = 0.25의 값을 사용했다.

식 (22)와 식 (23)의 가속도 벡터 u¨g와 속도 벡터 u˙g를 식 (21)에 대입하면, 다음과 같이 관성 효과가 반영된 강성행렬과 하중 벡터의 식으로 나타낼 수 있다.

식 (24a)로 부터 다음과 같이 변위-하중 관계를 정의할 수 있다.

최종적으로, 식 (20)의 변형률-변위 관계 행렬 B와 식 (25)의 변위-하중 관계 행렬 K^-1로 부터 다음 식 (26)을 도출할 수 있으며

여기서 행렬 A^가 변형률과 하중 사이에 직접 연관성을 갖는 관계 행렬이라 정의할 수 있다. 또한 행렬 A^는 최적의 센서 위치를 결정하기 위해 다시 사용되며 다음 절에서 이를 설명한다.

3.1장에서 계산한 행렬 A^를 D-최적 설계를 이용하여 공분산 행렬을 최소화시켜야 한다. 이를 위해 행렬 A^로부터 정의되는 정보 행렬 M^을 다음과 같이 계산하고,

Determinant 값이 최대화가 되는 행렬을 구축하기 위해 행 교환 과정을 반복하며 아래 식 (28)과 같이 행렬식을 반복적으로 계산한다.

식 (29)는 행이 추가된 A^ 정보 행렬의 행렬식을 도출하기 위해 새로운 후보 열 y^T이 추가되었을 경우 정보 행렬의 역행렬을 도출한다.

최종적으로 2장에서와 마찬가지로 행 교환 알고리즘의 반복을 통해 D-최적 설계를 얻음으로써 최소한의 오차로 하중 추정을 할 수 있는 변형 게이지가 부착되어야 하는 최적의 위치를 결정할 수 있다. 마지막으로, 변형 게이지로부터 얻은 입력 값들을 하중 벡터와 변형률 벡터의 관계식에 대입하여 추정하고자 하는 하중 벡터를 도출할 수 있다.

본 연구에서는 센서가 부착되지 않은 위치(non-sensor field)에서의 하중 추정을 하기 SEREP 기법을 응용하여 전 영역(full-field) 변위 및 하중 확장 행렬을 정의하였다. 확장 행렬은 소수의 센서가 부착되는 위치(sensor field)의 결과로 부터 전 영역의 하중과 변형률을 추정할 수 있다.

식 (26)의 변형률 벡터 ϵ_g는 구조물에 부착된 센서 위치에서 얻어진 변형률 벡터 ϵ_s와 센서가 부착되지 않은 위치의 변형률 벡터 ϵ_u로 구분할 수 있으며, 이를 다음과 같이 행렬 분할(matrix partitioning) 형태로 표할 수 있다.

여기서 Ψ_g와 Ψ_u는 변형률 벡터 ϵ_s와 ϵ_u에 대응되는 변형률 모드 기여 계수이다. 중요한 점은 변형률 벡터 ϵ_u는 센서 미부착 위치에서의 값으로, 실제로는 센서가 부착되지 않고 이론에 의해서 계산되는 변형률이다.

식 (30)로 부터 센서가 부착된 위치에 대한 변형률 벡터 ϵ_s는 다음과 같이 정의할 수 있다.

식 (31)의 양변에 ΨsT을 곱하면,

그리고 양변에 ΨsTΨs-1을 곱하면, 다음과 같이 정리될 수 있다.

최종적으로, 모드 기여 계수 벡터 q에 대한 식으로 정리를 한다면 식 (34)과 같이 표현할 수 있다.

구조물의 전체 영역(full-field)에서의 변형률 응답(strain response)를 도출하기 위해 식 (34)을 식 (30)에 대입하면, 다음 식 (35)와 같이 센서 위치에서의 변형률 계측값을 구조물 전체 자유도의 변형률 값으로 확장하는 변형률 확장 행렬(strain expansion matrix) E_s를 계산할 수 있다.

변형률과 마찬가지로 하중 벡터는 다음과 같이 행렬 분할(Matrix partitioning) 할 수 있다.

하중 추정 위치 자유도에 대한 하중 벡터는 식 (37)과 같이 정리할 수 있다.

마찬가지로 식 (37)도 식 (31)과 같이 미지수의 수가 방정식의 수와 같지 않으므로, 이 방정식의 경우 일반화 역행렬 방법(generalized inverse method)과 least-squares approach를 사용하여 풀어야 한다.

식 (38)을 정리하면 다음과 같이 표현된다.

식 (39)을 변형률 벡터 ϵ_g에 관한 식으로 정리를 하면

또한, 하중 벡터와 후보군 행렬도 변형률과 변형 모드 형상과 마찬가지로 일반화 역행렬을 사용하여 full-field force vector를 도출할 수 있는 Force expansion matrix E_f를 식 (42)과 같이 계산할 수 있다.

식 (42)에서 계산한 하중 벡터의 확장 행렬을 이용하여 통상적으로 하중이 가해지는 위치에 가해진 하중을 구조물 전체 노드 가해지는 하중으로 확장하여, 하중을 추정하고자 하는 위치를 사전에 설정하지 않고 구조물 전체 노드들에 가해지는 하중을 추정할 수 있다.

본 예제는 아래 Fig. 1과 같이 B = 37m, H = 87m, 두께 t = 0.025m인 자유도 155766개의 Jacket structure model을 Shell element로 모델링 하여 Main column 부분에 국부 종 방향 하중이 가해지는 동적 해석을 수행하였다.

Fig. 1 
The FE model of the jacket structure

유한요소 모델에 적용한 element type은 ABAQUS S4 shell element를 사용하였고, 재료는 탄성계수 E = 206000MPa, Poisson’s Ratio = 0.3인 연철(mild steel)을 사용하였다. 구조물에 가해진 동적 하중은 3초 동안 main column 한 개의 상단 끝단에 x축 방향으로 적용하였고, 가해진 하중 amplitude에 대한 load profile은 아래 Fig. 2와 같으며 jacket structure의 main column들의 밑단의 자유도들을 all-fixed 경계 조건을 적용하여 해석을 진행하였다.

Fig. 2 
Applied load profile of the jacket structure

Jacket structure FE 모델에서의 D 최적 설계를 통해 얻은 최적 센서 노드 위치를 아래 Fig. 3과 같이 설정하였고, 최적 센서 위치 노드와 센서 각도 방향은 Table. 1에 명시하였다.

Fig. 3 
Optimal gauge locations of jacket structure problem

Fig. 4는 target node에서 D-optimal design 기반 하중 추정법과 제안한 하중 추정법을 통해 추정한 전체 하중 벡터가 실제 FE 모델에 가해진 하중과 얼마나 일치하는지에 대해 도식하였다.

Fig. 4 
Applied load recovery results on target node (a) D-optimal inverse method (b) Proposed method

Target node에서의 하중 추정 오차율은 D-optimal 기반 하중 법은 24.7%이었고, 제안한 하중 추정법의 오차는 7.2%로 본 연구에서 제안한 하중 추정법의 정확도가 이전 기법의 정확도보다 우수한 것을 알 수 있다.

또한, 각 기법으로 추정한 하중을 이용하여 도출한 변위들과 실제 변위와의 일치 정확도를 Fig. 5와 같이 표현하였다. 각 하중 추정법을 이용하여 도출한 변위 값에 대한 오차율은 D-optimal 기반 하중 추정법은 34.3%이었고, 제안한 하중 추정법은 9.6%의 오차율로 변위에서도 하중과 마찬가지로 제안한 방법의 정확성을 검증하였다.

Fig. 5 
Displacement recovery results on target node (a) D-optimal inverse method (b) Proposed method

본 예제는 아래 Fig. 6과 같이 B = L = 1200mm, w = 400mm, l = 300mm 두께 t = 20mm인 자유도 110058개의 stiffened plate model을 shell element로 모델링 하여 bottom plate 부분에 y 방향 pressure 하중이 가해지는 동적 해석을 수행하였다.

Fig. 6 
The stiffened plate: (a) FE model (b) Section drawings of structural members

유한요소 모델에 적용한 element type은 앞선 예제와 마찬가지로 ABAQUS S4 shell element를 사용하였고, 재료는 탄성계수 E = 206000MPa, Poisson’s Ratio = 0.3인 연철(mild steel)을 사용하였다.

구조물에 가해진 동적 하중은 3초 동안 stiffened plate의 하부 Plate에 전체적으로 y축 방향 pressure 하중으로 적용하였고, 가해진 하중 amplitude에 대한 load profile은 아래 Fig. 7과 같으며 stiffened plate의 양 끝단 자유도들을 all-fixed 경계 조건을 적용하여 해석을 진행하였다.

Fig. 7 
Applied load profile of the stiffened plate

Stiffened plate 모델에서 진행한 D-최적 설계를 통해 얻은 최적 센서 노드 위치를 아래 Fig. 8과 같이 설정하였고, 최적 센서 위치 노드와 센서 각도 방향은 Table. 2에 명시하였다.

Fig. 8 
Optimal gauge locations of jacket structure problem

Fig. 9는 target node에서 D-최적 설계만을 이용하여 추정한 하중 벡터와 후보군 행렬의 모델 확장 기법을 통해 추정한 하중 값이 실제 ABAQUS 소프트웨어에서 FE 모델에 가해진 하중과의 일치 정확도를 나타내었다.

Fig. 9 
Applied load recovery results on target node (a) D-optimal inverse method (b) Proposed method

Target node에서의 하중 추정 오차율은 D-optimal 기반 하중 추정법은 18.4%이었고, 제안한 하중 추정법의 오차는 5.0%로 본 연구에서 제안한 확장 기법을 이용한 하중 추정법의 정확도가 이전의 기법의 정확도보다 향상된 것을 알 수 있다.

또한, 각 기법으로 추정한 하중을 이용하여 계산한 변위들과 실제 발생한 변위와의 일치 정확도를 Fig. 10과 같이 표현하였다. 각 하중 추정법을 이용하여 역 추정한 변위 값에 대한 오차율은 D-optimal 기반 추정법은 29.1%이었고, 제안한 하중 추정법은 13.6%로 변위에서도 하중과 마찬가지로 제안한 방법의 정확성을 검증하였다.

Fig. 10 
Displacement recovery results on target node (a) D-optimal inverse method (b) Proposed method