1. 서 론
해밀턴 이론(Hamilton, 1834; 1835)은 현재까지 동역학의 유일한 변분이론(variational principle)으로써, 가상일의 원리(principle of virtual work) 및 최소일의 법칙(principle of minimum total potential energy) 등의 다른 정역학의 변분이론과 같이 여러 분야에 널리 쓰이고 있다(Bretherton, 1970; Gossick, 1967; Landau, 1975; Slawinski, 2003; Tiersten, 1967). 하지만, 해밀턴 이론은 정역학의 변분이론들과는 다르게 두 가지의 큰 이론적 약점이 있다. 곧, 첫째, 해밀턴 이론은 동적시스템의 초기치(initial conditions)를 제대로 반영할 수 없고, 둘째, 에너지 비보존계(non-conservative systems)를 다룰 시 별도의 범함수(functional) Rayleigh's dissipation(Rayleigh, 1877)이 필요하다.
역사적으로 이런 해밀턴 이론의 단점을 극복하고자 한 대표적인 연구를 살펴보면, 크게 1) 합성적분(convolution)의 적용을 통한 Gurtin(1964a; 1964b)과 Tonti(1973)의 연구와 2) 정수차수가 아닌 미분(fractional derivative)을 적용한 연구(Riewe, 1996; 1997; Agrawal, 2001; 2002; 2008; Atanckovic et al., 2008; Baleanu, 2005; Cresson, 2007) 등이 있다. 하지만, 합성적분을 적용한 Tonti의 접근은 초기값이 0이 아닌 경우는 제대로 적용할 수 없는 어려움이 있고, Gurtin의 방법은 초기치 문제(initial-value problem)를 등가의 경계치 문제(boundary-value problem)로 환원시킴으로 완전한 강형의 미분방정식(strong form)을 회복할 수 없는 단점이 있다. 또한, 정수차수가 아닌 미분을 적용한 그 어떤 연구도 가장 간단한 경우인 선형탄성 단자 유도계에 관해서도 그 적용성이 검증된 바가 없다.
근래, 해밀턴 이론의 단점 모두를 극복한 혼합 합성 변분이론(Dargush and Kim, 2012)은 앞선 연구를 바탕으로 혼합 변수의 사용, 합성 적분의 사용, 그리고 정수차수가 아닌 미적분(fractional calculus)을 바탕으로 나온 것으로, 하나의 범함수인 action을 통해 에너지 보존 시스템과 에너지 비보존 시스템의 운동방정식과 초기조건을 이끌어 낸다.
곧, 혼합 합성 변분이론은 현재 동적 수치해석에 사용되는 반이산화방법(semi-discretization scheme-공간에 대해서는 일반적인 유한요소법 적용을 하고, 시간에 대해서는 여러 time integration techniques-Newmark's family, Finite difference methods, Runge-Kutta ethods, etc.-을 적용하는 방법)을 대신해 시공간(space-time)의 유한요소해석법을 직접 적용시킬 수 있는 이론적인 배경을 제공한다.
유의할 점은 혼합 합성 변분이론이 수학적으로 모델링되는 여러 가지 에너지 비보존 시스템을 아우르기에는 아직 그 일반적이고 구체적인 프레임이 성립되어 있지 않은 점이며, 이를 극복하는 연구가 현재 진행 중에 있다.
본 논문은 혼합 합성 변분이론을 활용한 시간에 대해 이차항의 형상함수가 적용된 세 가지 시간유한요소법을 제안하였고, 기본적인 수치해석법의 특징과 그 방법의 적절성을 주요 수치해석 예를 통해 검증하였다.
2. 본 론
이 장에서는 Fig. 1에서 보이는 선형탄성 단자유도계에 대한 혼합 합성 변분이론의 적용을 소개한 후, 이에 대한 수치해석법 개발 방법을 소개하고 이를 주요 예제를 통해 검증하고자 한다.
2.1 혼합 합성 변분이론의 선형탄성 단자유도계 적용
Fig. 1과 같은 선형탄성 단자유도계에서의 운동방정식은 다음의 식으로 정의되며,
(1)
초기조건은 다음과 같다.
(2)
식 (1)의 m, c, k, u, 
는 각각 질량, 감쇠 상수 (damping coefficient), 스프링 상수(stiffness), 질량의 변위, 외력을 의미한다. 또한, a 는 스프링 상수의 역수로 스프링의 연성(flexibility)을 나타내며, 위첨자 ∙ 는 시간에 대한 미분을 나타낸다.
혼합 합성 변분이론은 Fig. 1의 모델에 대해, 질량의 변위 u와 스프링내 힘의 충격량(impulse) J를 통해 다음의 범함수 action 
를 정의한다.
(3)
식 (3)에서 ∗은 시간간격 (0, t)에 대한 다음의 합성 적분(convolution)을 나타낸다.
(4)
또한, 식 (3)에서, 위첨자 ~ 은 다음의 left Riemann- Liouville semi-derivative를 나타낸다.
(5)
식 (5)의 Γ는 감마 함수(Gamma function)을 나타내며, Left Riemann-Liouville semi-derivative는 다음의 부분적분식(integration by parts)을 갖는다.
(6)
마지막으로, 식 (3)의 
는 초기의 하중 에 대한 충격량(impulse)으로 다음을 나타낸다.
(7)
곧, 는 초기시간 (t=0)으로 고려된 시간 이전(
로 표현됨)에 알려진 하중의 충격량을 의미한다.
혼합 합성 변분이론은 범함수 에 대한 정상성(stationarity)을 통해 식 (1), (2)를 완전한 약형(weak form)의 식으로 이끌어낸다.
먼저 범함수 에 대한 정상성은 다음의 식으로 표현된다.
δ=0; (8)
식 (8)의 첫 항 
을 비롯한 기존의 부분적분식이 적용가능한 부분은 다음과 같이 전개된다.
(9)
또한, ()을 비롯한 semi-derivatives을 포함하는 부분적분은 식 (6)을 따라 다음과 같이 전개된다.
(10)
식 (7)을 포함하는 모든 부분적분식을 전개하면, 다음의 식을 갖게 된다.
δ=0; (11)
시간간격 (0, t)내에 임의의 변분 δu 와 δJ를 고려하면, 식 (11)로부터 다음의 미분 방정식을 얻을 수 있다.
(12)
또한, 시간 t에서의 임의의 변분 
와 
를 고려하면, 다음의 초기조건이 성립함을 알 수 있다.
(13)
식 (11)에서 초기 변분 
은 주어진 초기값인 u0로 없어지며, 초기 변분 
는 주어진 초기값인 J0로 없어진다. 곧, 초기값인 u0, 
와 J0를 통해 초기값 J0는 식 (13)에 의해 산정되기에 이의 변분은 존재치 않게 된다.
식 (12)에서, 스프링 내의 힘 
이 다음의 식으로
(14)
표현됨을 감안하면, 혼합 합성 변분이론은 식 (1)은 물론, 재료 구성식(rate-constitutive relation)을 동시에 이끌어낸다. 덧붙여, 식 (13)에서 보듯, 혼합 합성 변분이론에서 초기조건은 독립적인 초기조건 
와 비독립적인 초기조건 
로 구성됨을 알 수 있다. 곧, 혼합 합성 변분이론은 하나의 범함수 에 대한 정상성을 통해, 강형의 미분 방정식과 초기조건 모두를 약형의 식으로 이끌어냄을 확인할 수 있다.
2.2 수치해석법 전개
식 (11)과 동등한 식 (8)은 시간에 대한 다양한 유한요소 해석법 개발에 근거를 제공한다. 곧, 변분항(variational term)을 가상 필드(virtual field)로 변분항이 아닌 항을 실제 필드(real field)로 간주하면, 식 (11)은 가상일의 원리가 혼합변수를 통해 동역학에 적용됨을 확인할 수 있다. 더불어, 식 (8)은 실제 필드 (δ, J)와 가상 필드 (δu, δJ) 모두가 균등하게 시간에 대한 C0 연속조건을 요구함을 확인할 수 있다.
Table 1
Developed algorithms
식 (8)을 통해 개발될 수 있는 여러 가지 동역학에 대한 해석법 중, 본 논문에서는 연속성 요구조건을 충분히 만족하는 시간에 대한 이차의 형상함수가 적용된 Table 1의 세 가지 시간유한요소해석법이 개발되었다.
2.2.1 시간에 대한 형상함수와 합성적분
시간간격(0
)에서 식 (8)의 주요변수는 시간에 대한 형상함수(shape functions)로 근사시킬 수 있다. 곧, 시간간격의 중심 
를 고려하면, 선형 형상함수 L과 이차 형상함수 Q는 다음과 같다.
(15)
(16)
(17)
(18)
(19)
또한, 식 (15)~(19)의 시간에 대한 미분은 다음과 같다.
(20)
(21)
(22)
(23)
(24)
식 (15)~(24)의 형상함수들의 시간간격(0)에서 합성적분(convolution)을 살펴보면 아래와 같다.
(25)
(26)
(27)
식 (25)에서 
의 합성적분이 어떻게 산정되었는지를 살펴보면, 다음과 같다.
(28)
형상함수의 Left Riemann-Liouville semi-derivative간의 합성적분은 아래의 멱함수(power functions)에 대한 Left Riemann-Liouville semi-derivative의 합성 적분식을 통해 산정될 수 있다.
(29)
곧, Left Riemann-Liouville semi-derivative간의 합성적분은 다음과 같다.
(30)
(31)
식 (29)을 활용한 식 (30)의 
계산 과정을 살펴보면 다음과 같다.
(32)
2.2.2 알고리즘의 전개방식
앞의 형상함수와 그에 따른 합성적분의 직접 산정을 통해, Table 1의 모든 알고리즘이 개발될 수 있는데, 이중 대표적으로 Uquad 알고리즘 개발방법을 소개하면 다음과 같다. 식 (8)의 
에서 
은 시간간격 
에서 다음의 매트릭스 형태로 표현된다.
(33)
식 (33)에서, 가상 필드와 실제 필드에 대한 아래첨자는 그 시간에서의 각 필드의 값을(discrete value) 나타내며, 합성적분 산정에는 합성적분의 특성 a * B = b * a이 고려되었다.
마찬가지로, 식 (8)의 다른 합성적분과 
는 아래의 식으로 표현된다.
(34)
(35)
(36)
(37)
(38)
(39)
식 (33)~(39)를 통해, 시간간격()에서 식 (8)은 다음의 식으로 표현된다.
δ=0 (40)
식 (40)에서 초기변분 δu0와 δJ0는 존재치 않는다: 주어진 초기조건 
으로부터 비독립적인 초기조건 
은 식 (13)으로 이끌어낼 수 있고, 이는 주요변수의 초기변분은 임의로 변할 수 없는 δu00, δJ00를 의미하게 된다. 곧, 식 (34)는 아래의 식으로 표현된다.
δ=0 (41)
식 (41)을 각 변분에(혹은, 가상 필드) 대해 정리하고 이 변분이 임의로 변하는 것을 주지하면, δ=0을 만족시키기 위해서는 다음의 식이 성립해야 한다.
(42)
(43)
(44)
식 (44)에서, uc는 다른 변수인 (J0, J1, u0, u1)로 표현될 수 있고, 이를 앞의 식 (42), (43)에 대입하면 변수 (J0, J1, u0, u1)로 표현되는 다음의 matrix equation을 얻을 수 있다.
(45)
식 (45)에서, X, Y, 그리고 Qu1은 다음을 나타낸다.
(46)
(47)
(48)
일반적인 시간간격 
곧, tn = nh에서, Uquad 알고리즘은 다음의 식으로 표현된다.
(49)
식 (49)에서, Qun은 다음을 나타낸다.
(50)
이와 마찬가지 방법으로 Jquad와 UJquad 알고리즘을 matrix equation으로 표현하면 아래와 같다.
Uquad algorithm: (51)
UJquad algorithm: (52)
식 (51), (52)에서, QJn과 QUJn은 다음을 나타낸다.
(53)
(54)
2.3 제안된 수치해석법의 특징
앞의 제안된 수치해석법은 다음의 간단한 matrix equation으로 나타낼 수 있다.
(55)
또한, 식 (51)은 동등한 아래의 식으로 표현된다.
(56a)
(56b)
식 (55)~(56)에서 매트릭스의 위첨자 -1은 역행렬을 의미한다.
외력과 감쇠가 존재하지 않는 에너지 보존 시스템(harmonic oscillator)인 경우, 식 (55)는 아래와 같다.
(57)
개발된 해석법의 안정성과 수치감쇠 발생여부는 식 (57)의 특이해(eigenvalues)의 크기에 의해 결정되는데, 각 알고리즘의 특이해는 Table 2와 같다.
Table 2에서 보듯, 모든 알고리즘에서 특이해의 크기는 다음의 식으로 표현된다.
(58)
곧, 제안된 알고리즘 모두는 time-step에 관계없이 무조건적인 안정성을 갖고, 수치 감쇠는 발생치 않음을 알 수 있다.
2.4 수치해석법 예제
제안된 수치해석법은 1) 감쇠 구조물의 자유진동(damped free vibration) 2) 조화하중에 의한 감쇠구조물의 진동(harmonic vibration of the damped system)의 두 가지 주요 예제릍 통해 검증되었다.
해석상 편의를 위해, 각 경우에 대해 질량 m=1, 스프링 강성 k=4π2, 감쇠 상수 0.2π에 3가지의 time-step cases(h=0.1, h=0.05, h=0.01)을 공통으로 하였고, 그 외 해석에 대한 주요정보는 Table 3과 같다.
Table 3
Numerical examples
| Damped free vibration | Forced damped vibration |
|---|---|
| u(0)=1
| u(0)=1
|
(0)=0
=0
Table 3에서처럼, 감쇠상수 c는 감쇠율(damping ratio) ε가 5%인 경우가 고려되었고, 각 해석결과는 정확한 해와 비교되었다. 제안된 해석법의 수렴특성을 확인하기 위해 추가로 OpenSees(Mckeena, F. et al., 2013)를 활용한 Newmark's methods의 해석도 수행하였다.
먼저, Newmark's constant average acceleration method에 의한 감쇠 시스템의 자유진동과 조화하중에 의한 결과는 Fig. 2~Fig. 3와 같다.
반면, 제안된 해석법에 의한 감쇠 시스템의 자유진동과 조화하중에 의한 결과는 Fig. 4~Fig. 9과 같다.
Fig. 2~Fig. 9을 비교시, 제안된 해석법 모두는 Newmark's constant average acceleration method에 비해 그 수렴의 속도가 현저히 빠른 것을 알 수 있다.
해석결과 비교에 사용된 Newmark's constant average acceleration method는 수치해가 time-step에 대해 무조건 안정적이기 때문에 Newmark's linear acceleration method에 비해 더욱 유용히 쓰이고 있다. 곧, Newmark's linear acceleration method는 다음의 조건을 만족시켜야하는 time-step에 대한 조건부 안정적인 해석법이다.
(59)
식 (59)의 Tn은 시스템의 고유주기를 나타내며, 단자 유도계가 아닌 경우, Newmark's linear acceleration method는 가장 작은 고유주기에 의해 수치해의 안정성을 만족시키는 time-step을 정해야 하는 단점이 있다. 수치해석 예제에 대한 Newmark's linear acceleration method를 통한 해석결과는 Fig. 10~Fig. 11과 같다.
Fig. 4~Fig. 11에서 보듯, 제안된 해석법은 Newmark's linear acceleration method보다 정해에 더욱 빠르게 수렴하는 것을 확인할 수 있다.
3. 결 론
본 연구에서는 동역학의 새로운 변분이론인 혼합 합성 변분이론에 근거해 시간에 대한 이차항의 형상함수가 적용된 해석법의 잠재력을 확인하고자 가장 간단한 선형탄성 단자 유도계로의 적용을 통해 그 해석법의 특징을 살펴보았다. 제안된 해석법은 time-step에 관계없이 항상 안정적으로, 에너지 보존 시스템의 경우 수치 감쇠가 발생되지 않는 에너지와 모멘텀이 보존되는 symplecticity property를 가지고 있음을 특이해 분석으로 확인할 수 있었다. 반면, 감쇠 시스템인 경우, Newmark's methods에 비해, time-step이 점점 작아질수록 정확한 해에 더욱 빠르게 수렴하는 것을 확인하였다. 곧, 본 연구를 통해 혼합 합성 변분이론을 토대로한 시간에 대한 고차형상함수를 활용한 해석법이 다자 유도계, 연속체를 비롯한 다양한 동적시스템 해석에 활용됨에 있어, 그 가능성과 잠재력이 확인되었다고 볼 수 있다.
제안된 해석법은 혼합 변수에 기반한 것으로 다자유도계와 연속체에 적용시, 기존의 변위기반의 해석법에 비해 계산비용 측면에 불리한 측면이 있지만, 혼합 변수에 기반한 해석법은 incompressible media의 동적 해석과 thermodynamics를 비롯해 다양한 적용이 가능한 장점 또한 갖고 있다. 향후, 이 기본 연구를 바탕으로 한 다양한 시공간 고차 유한요소해석법 개발과 그 적용을 기대해 본다.
