1. 서 론
구조물에 발생하는 균열은 구조물의 안전을 위협하는 주요한 원인 중의 하나로 다양한 구조물의 안전에 중요하게 고려해야 하는 문제이다. 그러나 균열문제는 변위의 불연속성이 나타나는 본질적인 어려움으로 인하여 일반적인 방법론으로는 해석하기가 쉽지 않다. 현재 많이 사용하고 있는 유한요소법은 균열이 진전할 때 요소망을 새로이 구성해야 하는 특별한 기법이 요구된다. 그러나 페리다이나믹스 이론은 지배방정식에서부터 변위에 대한 공간 미분항이 없고 적분방정식의 형태이므로 변위의 불연속성이 존재하더라도 균열의 진전을 자연스럽게 표현할 수 있다는 장점이 있다(Silling, 2000). 따라서 본 논문에서는 페리다이나믹 이론을 기반으로 균열문제에 대한 형상 설계민감도 해석법을 개발하였다.
이전의 많은 연구에서 동적 균열진전 문제에 대한 수치적 해석이 수행되었으나 설계민감도 해석은 유한요소법에 기반하여 정적인 문제에 대한 연구만 있었을 뿐 동적 문제에 대한 연구는 진행된 바 없었으며 균열진전 문제에 대한 설계민감도 연구는 균열의 제어와 최적설계라는 관점에서 중요하다. 실험적으로는 균열의 선단에 발생되는 응력을 제어하여 빠른 파괴단면을 얻으려는 연구가 있었고(Bowden et al., 1967) 최근에는 실험환경을 제어하면서 미소 스케일에서의 균열 패턴을 다양하게 만드는 연구(Nam et al., 2012)가 진행된 바 있다. 본 연구에서는 페리다이나믹스를 이용하여 실험적인 방법이 아닌 해석적 방법으로 형상변화에 대한 균열진전 문제의 설계민감도 해석을 수행하였으며 이는 페리다이나믹스에서 처음으로 시도된 연구이다.
설계민감도 해석법(Choi and Haug, 1983)은 지금까지 다양한 구조문제에 대해서 수학적으로 개발되어 왔는데, 변형의 정도가 큰 탄소성 재료의 다이나믹 문제에 대한 설계민감도 해석법(Cho and Choi, 2000)도 개발되었다. 특히 형상 설계민감도는 연속체에서 전미분(Material Derivative) 개념을 적용하여 물리적 영역 경계면의 수직방향으로의 변화만을 고려하여 얻어진다. 페리다이나믹스 이론에서는 형상 변화 중에 격자와 실제 형상이 함께 움직이므로 라그랑지안 관점으로 접근한다. 취성 탄성재료를 사용한 모델에 대하여 형상 설계민감도해석을 수행하고 그 정확성을 유한차분법(FDM)과 비교하였을때 매우 정확하고 효율적으로 얻을 수 있음을 확인하였다
2. 페리다이나믹스 이론
다음과 같이 영역 R에 속해있는 임의의 점 
에 대하여 페리다이나믹 지배방정식은 다음과 같이 구성된다.
(1)
여기서, Horizon H는 위치 벡터 x의 범위 내에 속하는 
를 포함하는 모든 주변 점들을 의미한다. 또한 u, b, 그리고 ƒ는 각각 변위 벡터, 체적력, 그리고 짝힘을 의미한다.
두 입자사이의 상대 기준위치벡터는 
로, 상대변위는 
로 정의한다. 따라서 변형상태에서의 상대 위치벡터는 
와 같다. 페리다이나믹스에서는 유한요소법과 같은 국부적인(local) 힘이 아닌 비국부적인(non- local) 힘이 작용하여 절점 x에 대한 일정한 범위 Hx 안에 들어오는 모든 주변 입자와의 짝힘을 고려한다. 본 연구에서는 Fig. 1과 같이 반경 의 원형 영역으로 그 범위를 가정하였다.
(2)
탄성재료에 대해서 포텐셜 ω를 이용하여 짝힘 ƒ는 다음과 같이 얻어진다.
(3)
여기서, ω는 한 결합 사이의 포텐셜을 의미하며 임의의 절점 x에서의 변형에너지 W는 다음과 같다.
(4)
특히 미소 탄성 취성재료(Ha and Bobaru, 2011)에 대하여 짝힘을 선형화하면 다음과 같이 표현된다.
(5)
여기서, C0은 선형 탄성계수로서 다음의 값을 이용한다.
(6)
E, v는 각각 영 탄성계수, 푸아송 비를 뜻한다. µ는 0 또는 1을 가리키는 함수로서 결합의 유무를 나타낸다. 결합이 끊어진 상태에서는 0이 되고 다시 1로 되돌아가지 못한다.
(7)
여기서, G0는 파괴에너지로서 단위 길이의 균열로 인하여 결합이 완전히 끊어지는데 필요한 에너지로 정의되며 실험적으로 얻어진다.
3. 페리다이나믹스의 이산화
페리다이나믹스에서 실제 해석영역은 작은 부피를 가지는 입자들로 이산화된다. 시간 스텝 n에서 임의의 i번째 입자에 대한 페리다이나믹 지배방정식 식 (1)의 이산화 형태는 관련된 부피 Vi와 결합된 주변 입자 j를 이용하여 다음과 같이 정리된다.
(8)
i 입자의 주변에서 짝힘을 가하는 입자들은 식 (2)를 만족하는 상대거리가 δ 인 원형 반경 안에 들어오는 입자들이다. 그러나 원형 반경의 경계 근처에 있는 입자들은 부피를 가지고 있기 때문에 실제 짝힘으로 계산되는 입자의 부피는 j 입자의 전체 부피보다 작아야 한다. 이에 따라 체적율(Volume Fraction) vi이 도입된다.
(9)
실제 페리다이나믹에서의 짝힘을 계산하기 위해서는 원형반경 안에 들어오는 짝힘에 영향을 주는 j 입자의 부피를 계산해야 한다.
예를 들어 p입자는 j입자의 원형 반경 안에 완전히 들어오지만 q입자는 j입자와 일부 부피만 영향을 주게 되는데 이에 대한 고려는 페리다이나믹스에서는 수치적 근사로 해결하고 있다. 실제로 다양한 알고리즘이 쓰이고 있는데 본 연구에서 정사각형으로 격자가 이루어지는 경우에 대하여 아래(Hu Wenke, 2012)의 알고리즘을 활용하였다.
(10)
4. 형상 설계민감도 해석
Fig. 3와 같은 영역의 변화를 살펴보자. 변환 관계는 
, 
로 주어지고 τ는 설계 변환 파라미터를 뜻한다.
T는 테일러 급수를 이용하여 선형화된 관계를 생각한다.
(11)
여기서, 속도장 λ는 다음과 같이 정의된다.
(12)
이를 이용하여 지배방정식 식 (1)에 전미분을 취하면 다음과 같이 정리된다.
(13)
여기서, 
를 의미한다. 우변의 첫 번째 항은 짝힘에 대한 전미분 항이다. 이는 η와 ε에 대한 식으로 분해된다. 또한 Horizon 경계 부근의 체적율을 고려하기 위해서 Vs를 이용하여 식 (13)의 우변 첫 번째 항을 다시 정리하면 다음과 같다.
(14)
식 (14)에서 첫 번째 항은 짝힘에 관한, 두 번째 항은 체적율에 대한 항이다. 식 (14)에 대한 이산화 식은 아래와 같다.
(15)
여기서, 
, 
는 각각 상대 기준위치, 상대변위에 대한 미분항으로서 다음의 성질을 만족한다.
특히 선형화된 짝힘 식 (5)를 이용하면 
, 
는 다음과 같이 정리된다.
(16)
(17)
5. 애조인 변수법
, 
에 대해서 식 (1)에 애조인 변수 λ(x,t)를 곱하면 다음과 같이 정리된다.
(18)
식 (18)을 τ에 대해서 변분을 취하면 다음과 같다.
(19)
부분적분을 이용하여 식 (19)에 다음의 관계를 이용한다.
(20)
여기서, 형상변화의 초기조건은 0이다.
목적함수를 마지막 시간에서의 영역 적분값과 전체 시간에 대한 적분 값으로 일반화 하면 다음과 같이 쓸 수 있다.
(21)
식 (21)을 부분적분을 이용하여 변분을 취하고, 식 (19)을 대입하면 λ에 대한 목적함수의 설계민감도를 얻을 수 있다.
(22)
λ를 얻기 위한 애조인 시스템은 다음과 같이 구성된다.
(23)
그리고 경계조건은 다음과 같다.
(24)
(25)
동적 문제에서 애조인 방정식은 종기치 문제(terminal value problem)이므로 초기조건은 시스템의 가역성을 이용하여 얻어진다. 가역적 시스템은 2계 미분 방정식으로 정의될 수 있고, 시간과 속도의 반대방향에 대해서 변하지 않는 성질을 갖는다. 일반적인 동적 시스템에 대해서 해석과 애조인 시스템은 시간 대칭성을 갖지만 설계민감도 방정식은 경로에 의존하는 힘이 존재하여 시간에 대한 대칭성을 가지지 못한다(Jang et al., 2014). 해석과 애조인 시스템이 시간에 대한 대칭성을 가진다는 것은 시간 반대방향으로 적분할 때 정확한 해석 경로를 따라갈 수 있고 해석의 경로를 저장할 필요가 없기에 매우 효율적으로 애조인 설계민감도를 계산할 수 있다. 그러나 페리다이나믹스에서는 결합 파라미터 µ로 인하여 해석에서의 시간 대칭성을 기대하기는 어렵다. 따라서 애조인 시스템을 풀기 위해서는 저장된 해석값이 불가피하게 필요하다.
6. 병렬 기법
큰 스케일의 모델링과 계산 비용의 효율성을 높이기 위해서 병렬 기법에서 중요한 것은 각각의 프로세서에 할당되는 부하의 균형(load balancing)과 통신을 최소화하는 것이며 본 연구에서는 이진분해(Binary decomposition; Bowden et al., 1967)을 이용하여 페리다이나믹 해석과 설계민감도에 적용하였다. 이진분해에서 영역을 각 프로세서에 할당하기 위한 분할법은 순차적으로 영역크기가 거의 비슷하도록 두 부분으로 나누는 것이다. Fig. 4와 5에서는 임의의 영역을 8개의 프로세서에 할당하도록 나누는 과정이다.
영역 분할 후에는 통신이 필요한 주변 영역의 프로세서를 결정해야 한다. 이를 위해 각 입자의 위치마다 Horizon 안에 들어오는 주변 입자를 모두 찾아서 주변 프로세서를 결정해준다. 이 과정은 O(NXN)(N : 모든 입자의 개수)의 계산비용을 필요로 하여 계산 비용이 높지만 해석과정 중에 처음 단 한번만 수행해 준다. 본 연구에서는 MPI 라이브러리를 활용하였고, 클러스터 CPU는 Intel Xeon E5-2600로 16개의 core를 이용하였다.
7. 수치 예제
Fig. 6와 같이 a=0.01m 길이의 초기 균열을 가지는 직사각형 평판을 고려하자. 재료의 물성치는 다음과 같다; E=720×108(Pa), ρ= 2440(kg/m3), G0=135. 본 예제 모델에서는 푸아송 비는 v=1/3으로 설정하였다. 또한 크기는 0.025m×0.01m이고 위쪽과 아래쪽 방향으로 작용하는 균일한 힘(3MPa)을 가하였다. Velocity-Verlet 알고리즘을 사용하였고, 시간 간격은 0.01µs이다. 전체 영역은 144,000개의 이산화된 격자로 나누었다.
목적함수는 마지막 시간에서의 특정한 위치의 변위로 선택하였고, 설계변수는 하중조건이 가해지는 면은 바뀌지 않는 범위에서 오른쪽 끝이 볼록하게 변화하도록 설정하였다.
(26)
먼저 애조인 변수법(AVM)과 직접미분법(DDM)으로 계산한 형상 설계민감도 값의 시간 비교를 한 결과를 Table 1에서 보여준다. 1개의 설계변수에 대한 민감도 해석을 수행하였을 때는 애조인 시스템을 풀어야하는 문제때문에 AVM이 DDM보다 시간이 더 걸리지만 설계변수가 많을수록 DDM에 비해서 AVM이 훨씬 효율적인 방법임을 보여준다. 이는 설계변수가 많아도 애조인 시스템은 단 한번만 수행하면 여러 설계변수에 대한 설계민감도를 얻을 수 있기 때문이다.
Table 1
Comparison of computation(AVM DDM)
| Analysis | 1 Design variable | 4 Design variables | |||
|---|---|---|---|---|---|
| AVM | DDM | AVM | DDM | ||
| Time (sec.) | 740.566 | 1057.467 | 778.881 | 1090.117 | 3215.443 |
| Normalized | 1.000 | 1.428 | 1.052 | 1.472 | 4.342 |
형상 설계민감도 해석에서는 체적율(volume fraction)에 대한 미분이 필요하므로 적당한 체적율의 선택이 중요하다. 기존의 연구에서 일반적으로 쓰이는 부분적 선형함수(piece- wise linear function)의 volume fraction, 식 (10)은 형상 설계민감도 해석에 적당하지 못하다. 식 (27)은 모든 구간에서 미분 가능한 체적율 함수이다.
(27)
시뮬레이션은 식 (10)을 쓴 결과(Fig. 7)와 식 (27)을 쓴 결과(Fig. 8)가 거의 차이나지 않는다.
그러나 민감도 해석결과는 비슷하지 않다. 두 결과 모두 대부분의 위치에서 유한차분법(FDM)과 거의 일치하는 결과를 보이지만 일부 위치(x-disp(A))에서의 목적함수에 대한 설계민감도 값이 다르다. 이는 유한 차분법보다는 해석적 설계민감도에서 생기는 체적율(volume fraction)에 의한 오차라고 추정된다. 식 (27)을 적용했을 때 설계민감도는 유한차분법과 훨씬 잘 일치하는 결과를 보인다.
Table 2
Comparison of design sensitivity using piece-wise volume fraction
Table 3
Comparison of design sensitivity using quadratic volume fraction
2차 체적율(quadratic volume fraction) 함수를 이용하여 애조인 변수법으로 유도한 형상 설계민감도 값을 유한차분법(FDM)과 비교하여 정확성을 확인하였다. 목적함수는 식 (26)과 같고 이를 이용하여 유도되는 목적함수에 대한 설계민감도식은 다음과 같다.
(28)
그리고 애조인 시스템과 그에 대한 경계조건은 다음과 같이 유도된다.
(29)
(30)
λ(tT)=0 (31)
Table 4
Comparison of adjoint design sensitivity
Table 4에서는 애조인 변수법으로 계산한 설계민감도를 유한차분법과 비교하였다. 목적함수의 A, B, C 위치는 Fig. 7, Fig. 8에서의 위치와 동일하다. 유한차분법에서의 설계변수의 변화량은 원래 형상의 0.001%로 정하였다. 유도된 형상설계민감도 해석값을 유한차분과 비교하였을 때 매우 정확하게 얻어짐을 확인하였다.
8. 결 론
본 연구에서는 페리다이나믹스 이론을 이용해서 균열진전문제에 대한 형상 설계민감도 해석을 수행하였다. 페리다이나믹스 이론은 균열 때문에 야기되는 변위의 불연속성을 표현하는데 특별한 기법 없이도 자연스러운 전개가 가능하고 초기 균열도 유한요소법보다 정확한 모델링이 가능하다. 설계민감도는 애조인 변수법을 이용하여 전개하였기 때문에 많은 설계변수가 존재하더라도 매우 효율적인 설계민감도 해석이 가능하다. 그러나 애조인 시스템을 풀기 위해서는 페리다이나믹 모델에 따라 해석경로를 저장해야 한다. 형상 설계민감도 해석에서 체적율의 영향을 고려해야 하는데 이에 따라 C1-연속성을 만족할 때 더 정확한 설계민감도를 얻을 수 있었다. 페리다이나믹스와 설계민감도 해석에 병렬기법을 적용하여 효율적이고, 큰 스케일의 모델 해석이 가능하였다.
