Journal of the Computational Structural Engineering Institute of Korea. 2015. 425-431
https://doi.org/10.7734/COSEIK.2015.28.4.425

ABSTRACT


MAIN

1. 서 론

파괴역학은 충격, 균열 진전문제 등의 공학적 분야에서 다양하게 이용된다. 지난 수십 년간 실험적으로 동적 하중에 대한 재료의 손상 또는 파괴에 대한 연구가 진행되어 왔다. 실험적으로는 유리같이 균열의 제어가 어려운 재료를 이용하여 균열이 분기되지 않고 균열을 빠르고 깔끔한 단면으로 잘리도록 하는 외부 힘을 제어하려는 연구가 진행된 적이 있고(Bowden, et al., 1967), 최근에는 다양한 외부환경을 변화시켜 보면서 미소 균열의 패턴을 변화시킨 연구도 진행되었다(Nam, et al., 2012). 그러나 본 연구에서는 실험적 방법이 아닌 최적화 기법을 이용하여 균열의 제어 문제를 이론적으로 접근하였다.

이론적으로 균열문제는 해석영역 안에서 변위의 비연속이 발생하기 때문에 일반적인 연속체 이론을 적용하기 어렵다. 이러한 문제를 해결하기 위해 X-FEM 기법이 도입되는 등 연속체이론 기반에서 해석하려는 이전의 노력이 있었다(Belytschko and Black, 1999). 그러나 본 연구에서는 연속체 기반 해석이론이 아닌 페리다이나믹스 이론을 이용하여 접근하였다. 연속체 기반에서 특별한 방법들을 도입하여 비연속적인 해를 해석하는 다른 방법과 달리 페리다이나믹스에서는 자연스럽게 비연속해를 얻을 수 있다는 장점을 가지고 있다. 따라서 본 연구에서는 페리다이나믹스 이론을 균열문제의 해석기법으로 채택하였다.

설계민감도법은 재료, 유체, 열, 비선형 문제뿐만 아니라 변형정도가 큰 동역학 문제에 대해서도 개발되었다(Cho and Choi, 2000). 특히 애조인 변수법은 설계변수의 개수가 많을 때 계산비용을 매우 절감할 수 있는 방법으로서 최적화기법에 많이 이용된다. 본 연구에서는 재료의 밀도를 설계변수로 하여 재료의 변형에너지에 대한 설계민감도를 애조인 변수법으로 계산하였다. 이를 활용하여 간단한 최적설계 예제를 통해 균열의 진행을 이론적인 접근방법으로 제어가 가능하다는 것을 확인하였다.

2. 페리다이나믹스 이론

다음과 같이 영역 에 속해 있는 임의의 점 에 대하여 페리다이나믹 지배방정식은 다음과 같이 구성된다.

ρu..(x,t)=Hf(u~(x~,t))-u(x,t),x~-x)dVx~+b(x,t)      (1)

여기서, H는 위치 x의 범위 안에 들어오는 다른 입자의 위치 https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-04/11TK082015280401/images/10.7734.28.4.425.F152.png를 포함하는 모든 주변 점들을 의미한다. 또한 u, b, 그리고 f는 각각 변위 벡터, 체적력, 그리고 짝힘을 의미한다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-04/11TK082015280401/images/10.7734.28.4.425.F102.png
Figure 1

Influenced region and displacement at x

두 입자사이의 상대 기준 위치벡터를 https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-04/11TK082015280401/images/10.7734.28.4.425.F103.png로 정의하고, 상대변위도 마찬가지로 https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-04/11TK082015280401/images/10.7734.28.4.425.F104.png로 정의한다. 변형된 입자 사이의 상대 위치벡터는 ξ+η와 같다. 페리다이나믹스 이론에서는 유한요소법과 같은 국부적(local)힘이 아닌 비국부적(non-local)힘이 작용하여 절점 x에 대한 일정한 소영역(horizon) Hx 안에 들어오는 모든 주변 입자와의 짝힘을 고려한다. 본 연구에서는 Fig. 1과 같이 반경 의 원형 영역으로 그 범위를 가정한다.

Hx= { x R "  | x-x~ | }       (2)

포텐셜 를 이용하면 짝힘 는 다음과 같이 얻어진다.

f(η,ξ)=w(η,ξ)ηη,ξ      (3)

ω는 한 결합사이의 포텐셜을 의미하며 변형에너지 밀도로 이해할 수 있다. 임의의 절점 x에서의 변형에너지 W는 다음과 같다.

W=12Hxw(η,ξ)dVξ      (4)

특히 미소탄성 취성 재료(Ha and Bobaru, 2010)에 대하여 짝힘을 선형화하면 다음과 같이 표현된다.

f(η,ξ)=η+ξ | | η+ξ | | μCos      (5)

여기서, C0은 선형 탄성계수로서 다음의 값을 이용한다.

Co=6Eπδ3(1- υ )      (6)

E, ν는 각각 영 탄성계수, 푸아송 비를 뜻한다. μ는 0 또는 1을 가리키는 함수로서 결합의 유무를 알려준다. 결합이 끊어진 상태에서는 0이 되고 다시 1로 되돌아가지는 못한다.

s0=4πG09Eδ      (7)

G0는 파괴에너지로서 단위 길이만큼 균열로 인해 결합이 완전히 끊어지는데 필요한 에너지로 정의되며 이 값은 실험적으로 얻어진다.

3. 이산화

페리다이나믹스 이론에서 실제 해석영역은 작은 부피를 가지는 입자들로 이산화된다. n번째 시간 스탭에서 i번째 입자에 대한 페리다이나믹 지배방정식 식 (1)의 이산화 형태는 관련된 부피 Vi와 결합된 주변 입자 j를 이용하여 아래와 같이 정리된다.

ρu..inVi=jFif(ujn-uin ,   xj-xi)Vj~Vi+binVi      (8)

i입자의 주변에서 짝힘을 가하는 입자들은 (2)를 만족하는 상대거리가 δ인 원형 반경 안에 들어오는 입자들이다. 그러나 원형 반경의 경계 근처에 있는 입자들은 부피를 가지고 있기 때문에 실제 짝힘으로 계산되는 입자의 부피는 j입자의 전체 부피보다 작아야 한다. 이에 따라 부피율(volume fraction) νs이 도입된다.

Vj~= υ s(xj-xi)Vj      (9)

페리다이나믹스 이론에서의 짝힘을 계산하기 위해서는 소영역 안에 들어오는 짝힘에 영향을 주는 j입자의 부피를 계산해야 한다. 따라서 소영역의 크기는 페리다이나믹스 이론에서 해석결과에 많은 영향을 미치므로 중요하게 다루어진다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-04/11TK082015280401/images/10.7734.28.4.425.F114.png
Figure 2

Horizon of particle on discrete grids

예를 들어 p입자는 i입자의 원형 반경 안에 완전히 들어오지만 q입자는 i입자와 일부 부피만 영향을 주게 되는데 이에 대한 고려는 페리다이나믹스에서 수치적 근사로 해결하고 있다. 실제로 다양한 알고리즘이 쓰이고 있는데 본 연구에서 정사각형으로 격자가 이루어지는 경우에 대하여 다음의 식(Hu Wenke, 2012)을 활용하였다.

υ s(r)= { 1, | |  r  | |  δ -2-1 | | r | | +(δ+12),δ-2 | |  r  | | δ+20,δ+2 | |  r  | |       (10)

4. 애조인 변수법

페리다이나믹 지배방정식이 모든 시간 동안 라그랑지 멀티플라이어의 곱에서도 식이 만족한다고 하면 다음과 같다.

0tTΩλT [ ρz..-Hf(z~-z,x~-x)dVx~-b ] dVdt=0      (11)

그리고 애조인 변수 λ가 설계변수에 대해 독립이므로 전체 식에 대하여 전미분을 취하면 다음과 같이 쓸 수 있다.

0tTΩλT [ ρz..'-ρ,uδuz..-b,uδu ] dVdt-      (12)
0tTΩλT [ H(fz~z'~+fz~z'+ f,  uδu)dVx~ ] dVdt=0

식 (12)을 부분 적분함으로써 다음의 관계를 얻는다.

0tTλTρz'..dVxdt      (13)
= [ Ω(λTρz'.-λ.Tρz')Vx ] 0tr+0trΩz'Tρλ..dVxdt

여기서, fz에 대해서 다음과 같은 대칭성을 가진다.

fzj=-fzi      (14)

또한 이를 이용하여 다음의 교환관계가 성립함을 알 수 있다.

0trλT(Hfz~z'~Vx~)dVxdt      (15)
=0trz'T(Hfz~z'~Vx~)dVxdt

δuz'https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-04/11TK082015280401/images/10.7734.28.4.425.F152.png에 대하여 독립적이기 때문에 식 (13)은 다음과 같이 정리된다.

0trΩλT [ ρ,uz..-b,u-Hf,udVx~)dudVdt      (16)
+Ωρ(λTz'.-λ'Tz')dVx | tr+0trΩz'T [ ρλ..-H(fz~λ~+fzλ)dVx~ ] dVdt=0

z'https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-04/11TK082015280401/images/10.7734.28.4.425.F130.png는 초기조건에서 설계변수에 무관하므로 https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-04/11TK082015280401/images/10.7734.28.4.425.F130.png(x, 0;u)=0, z'(x, 0;u)=0. 주어진 종료시간(terminal time)에 대하여 과도동역학(Transient dynamics)에서의 일반적인 성능함수(Performance measure)의 형태는 다음과 같다.

ψ=Hg(z,z.)dVx | tr+0trHh(z,z.)dVxdt      (17)

여기에 설계변수에 관한 변분(variation)을 취하고 정리를 하면, 다음과 같은 형태를 얻는다.

ψ'=H [ g,uδu+(g,z+h,z)z'+g,z.z'. ] dx | tr      (18)
+0trH [ h,uδu+(h,z-(h,z))z' ] dxdt

식 (18)에서 앞서 구한 식 (16)을 빼면 애조인 변수법을 이용한 목적함수 ψ에 대한 설계민감도를 얻을 수 있다.

ψ=Ω [ g,uδu ] dx | tr      (19)
+0trH [ h,u-λ { ρ,uz..-b,u-Hf',udx~ } ] δudxdt

여기서, λ는 애조인 변수로서 이를 구하기 위한 애조인 시스템은 다음과 같이 구성된다.

ρλ..-H(fz~λ~+fzλ)dVx~=h,z-dh,z.dtT      (20)

그리고 종료조건은 다음과 같이 주어진다.

ρλ.(tT)=-(g,z+h,z.)T | tT      (21)
ρλ(tT)=g,z.T | tT      (22)

식 (20)~(22)의 애조인 시스템은 종료시간에서의 경계 조건을 가지는 종기치 문제(terminal value problem)이고 시간에 대해 역방향을 가지는 동적 시스템임을 의미한다. 일반적인 2계 미분 동역학 문제에서의 애조인 방정식은 원래의 시스템과 시간에 대한 대칭성을 가진다(Jang, et al., 2014). 그러나 페리다이나믹 시스템에서는 로 인해 대칭성을 잃게 되므로 애조인 시스템을 위해 해석에서의 값을 저장해 두어야 한다(Kim and Cho, 2014).

5. 최적설계 기법

본 연구에서는 구조물의 균열진전 방향을 제어하는 것을 목적으로 하며 균열을 해석하기 위해서는 구조물에 걸리는 응력, 또는 변형에너지를 계산하여야 한다. 페리다이나믹스 이론에서는 구조물 입자 사이의 거리가 임계값을 넘어가면 결합(bond)이 끊어지면서 균열이 발생한다. 그리고 결합이 끊어지는 임계거리는 변형에너지의 실험값과 비교하여 결정된다. 또한 균열의 분기는 응력의 집중에 따라 달라진다는 것이 이전 실험논문에 의해 알려져 있다. 본 논문에서는 균열이 발생할 때 원하는 위치의 변형에너지를 최소화시켜 균열의 분기를 조절하는 방향으로 최적설계를 진행하였다.

구조물의 변형에너지는 식 (4)와 같이 정의된다. 이를 이용해서 목적함수를 다음과 같이 정의할 수 있다.

       Minimize ψ= [ Ω-12ωdΩ ] T      (23)

       Subject to ΩρdVM      (24)

ρminρρmax      (25)

여기서, https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-04/11TK082015280401/images/10.7734.28.4.425.F141.png는 변형에너지를 줄이고자 하는 목적 영역을 나타내고, ω는 변형에너지 밀도를 의미한다. 선형 탄성재료에 대하여 ω는 다음과 같이 정의된다.

ω=H12cs2 | ξ | dV      (26)

제약조건으로는 전체영역의 질량을 일정량 이하로 설정하였으며 약 80개의 부분영역으로 나누어 매개화하였다.

6. 수치 예제

Fig. 3과 같은 직사각형 평판에서 a=0.01m의 길이를 가지는 초기 균열을 고려하자. 본 예제 모델에서는 푸아송 비는 ν=1/3으로 설정하였다. 또한 크기는 0.025m×0.0125m이고 위쪽과 아래쪽 방향으로 작용하는 균일한 힘(3MPa)을 가하였다. 사용된 재료의 물성치는 다음과 같다.

E=720X108(Pa) , ρ=2440(kg/m3), G0=135

Velocity-Verlet 알고리즘을 사용하였고, 시간 간격은 5ns이며 전체영역은 5,000개의 이산화된 격자로 나누었다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-04/11TK082015280401/images/10.7734.28.4.425.F143.png
Figure 3

Rectangular plate with an initial crack

소영역 크기(Horizon size) δ는 일정한 간격을 가지는 각 격자의 폭인 x의 3배로 설정하였다. 소영역 크기는 작아질수록 파괴되는 영역도 작아지기 때문에 균열의 경로에 매우 중요한 요소 중 하나이다. 본 연구에서는 이전 연구의 결과를 참조하여 결정하였다(Ha and Boraru, 2010).

Fig. 4Fig. 5에서는 해석이 진행되면서 3.1μs와 10μs에서의 각 입자의 파괴정도와 변형 에너지를 보여주고 있다. 그림에서 확인할 수 있듯이 균열선단 주변에 있는 입자는 매우 높은 변형에너지를 저장하고 있게 됨을 확인할 수 있다. 따라서 균열선단 주변 변형에너지의 양을 줄일 수 있다면 균열의 진전 방향을 조절할 수 있을 것이다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-04/11TK082015280401/images/10.7734.28.4.425.F144.pnghttps://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-04/11TK082015280401/images/10.7734.28.4.425.F145.png
Figure 4

Crack propagation(Top: 3.1μs, Bottom: 10μs)

https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-04/11TK082015280401/images/10.7734.28.4.425.F146.pnghttps://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-04/11TK082015280401/images/10.7734.28.4.425.F147.png
Figure 5

Strain energy(Top: 3.1μs, Bottom: 10μs)

그에 앞서 유도된 설계민감도 해석결과를 Table 1에서 보여준다. Table 1에서는 변위에 대해 애조인 변수법으로 계산한 설계민감도 값을 유한차분법과 비교하였고 Table 2에서는 변형에너지의 설계민감도를 유한차분법과 비교하였다. 애조인 변수법으로 계산한 설계민감도는 균열이 진전되더라도 유한차분법으로 계산된 값과 잘 일치함을 알 수 있다.

Table 1

Design sensitivity of displacement

time node AVM FDM Accuracy
620 700-y 1.20515E-09 1.20545E-09 99.975%
290-y -5.27507E-11-5.27561E-1199.990%
2000700-y -4.63696E-07-4.63696E-07100.000%
290-y 1.00230E-071.00230E-07100.000%
Table 2

Design sensitivity of strain energy

time node AVM FDM Accuracy
620 700 2.21076E+032.21077E+0399.999%
290 1.25825E+051.25832E+0599.995%

최적설계를 진행하기 위해서 Fig. 6과 같은 인장응력 하의 모델을 고려한다. 평판의 크기는 Fig. 3과 동일하고 A로 표시된 위치가 변형에너지를 최소로 하도록 하는 목적함수의 위치이다. 이는 Fig. 5에서 볼 수 있듯이 3.1μs에서 변형에너지가 가장 큰 지점이며 Fig. 4에서와 같이 분기가 발생되는 위치이다. 이 위치의 변형에너지를 줄이는 최적설계를 진행하여 균열의 분기 위치를 조정하여 균열진전을 제어하는 최적의 재료분포를 얻고자 한다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-04/11TK082015280401/images/10.7734.28.4.425.F148.png
Figure 6

Crack propagation model for optimal design

최적설계에 필요한 ρmin은 1952, ρmax는 3416 그리고 M은 0.5로 최적화 모델을 선택하였다. 또한 설계변수에 대한 매개화는 Fig. 7과 같이 80개의 격자로 모델링하였다. Table 3에서는 초기의 A위치에서의 변형에너지와 최적화가 끝났을 때의 변형에너지를 보여주고 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-04/11TK082015280401/images/10.7734.28.4.425.F149.png
Figure 7

Parameterization of design variables

Table 3

Comparison of strain energy

DesignInitialOptimal
Objective function20316.98179311.1687

최적 재료분포로 만들어지는 변형에너지는 초기 값에 비해서 약 45.8%로서 감소한 것을 확인할 수 있다. 구조물의 균열진전 형태는 균열선단에서의 응력의 크기에 따라 분기가 달라진다. 따라서 이를 확인하기 위해 최적 밀도분포로 같은 해석 모델을 이용해 10까지 해석을 수행하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-04/11TK082015280401/images/10.7734.28.4.425.F150.pnghttps://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-04/11TK082015280401/images/10.7734.28.4.425.F151.png
Figure 8

Crack propagation before(Top) and after(Bottom) Optimization

Fig. 8에서 확인할 수 있듯이 원래의 분기가 일어나는 A점보다 분기의 위치가 더 뒤로 이동한 것을 확인할 수 있다. 이것은 균열선단에서의 응력의 크기가 분기를 결정한다는 이전의 연구 결과와 같다(Bowden, et al., 1967).

7. 결 론

본 연구에서는 페리다이나믹스 이론을 이용해서 균열진전문제에 대한 형상 설계민감도 해석을 수행하고 이를 활용하여 최적설계를 수행하였다. 균열진전 해석을 위해서 적용한 페리다이나믹스는 비국부적(non-local) 시스템으로서 기존의 연속체 이론에서 불연속성을 표현하기 어려웠던 문제를 자연스럽게 해결할 수 있는 방법이다. 설계민감도 해석에서는 애조인 변수법을 이용하여 매우 효율적이고 정확한 설계민감도 해석을 수행하였다. 특히 균열이 진전되어 불연속성이 많아지더라도 유한차분법과 비교하여 정확한 설계민감도 값을 얻을 수 있음을 확인하였다. 다만 페리다이나믹에서는 애조인 시스템이 해석과 시간 대칭성을 가지지 못하기 때문에 애조인 시스템을 풀기 위해서는 페리다이나믹 해석에서의 계산과정을 저장하고 있어야 한다. 인장응력을 가지는 간단한 평판 구조물에 대해 분기되는 위치의 변형에너지를 최소로 만드는 재료밀도 분포 형태를 최적설계 기법으로 얻었다. 얻어진 재료밀도 분포를 이용하여 페리다이나믹 해석을 수행하였을 때 원래의 균열 분기점이 뒤로 많이 밀려나 분기의 위치를 제어할 수 있음을 확인하였다.

Acknowledgements

이 논문은 2013 년도 미래창조과학부의 재원으로 한국연구재단의 지원을 받아 수행된 연구입니다(No. 2010-0018282). 저자들은 연구비 지원에 깊은 감사를 드립니다.

References

1
T., Belytschko and T., Black, 1999. Elastic Crack Growth in Finite Elements with Minimal Remeshing. Int. J. Num. Methods Eng., 45, pp.601-620.
10.1002/(SICI)1097-0207(19990620)45:5<601::AID-NME598>3.0.CO;2-S
2
F.P., Bowden, J.H., Brunton, J.E., Field and A.D., Heyes, 1967. Controlled Fracture of Brittle Solids and Interruption of Electrical Current. Nature, 216, pp.38-42.
10.1038/216038a0
3
S., Cho and K.K., Choi, 2000. Design Sensitivity Analysis and Optimization of Non-linear Transient Dynamics Part I-sizing Design. Int. J. Numer. Methods Eng., 48, pp.351-373.
10.1002/(SICI)1097-0207(20000530)48:3<351::AID-NME878>3.0.CO;2-P
4
Silling S.A. 2000Reformulation of Elasticity Theory for Discontinuities and Long-range ForceJ Mech Phys Solid481752092810.1016/S0022-5096(99)00029-0
5
Ha Y.D Boraru F. 2010Studies of Dynamic Crack Propagation and Crack Branching with Peridynamics Int. J. Fracture162229244
10.1007/978-90-481-9760-6_18
6
H.-J., Jang, J.-H., Kim, Y., Park and S., Cho, 2014. Adjoint Design Sensitivity Analysis of Molecular Dynamics in Parallel Computing Environment,. Int. J. Mech. & Mater. Design,, 10, pp.379-394.
10.1007/s10999-014-9253-2
7
Hu W. 2012Peridynamic Models for Dynamic Brittle FractureEngineering Mechanics Dissertations & Theses,Paper 28
8
K.H., Nam, I.H., Park and S.H., Ko, 2012. Patterning by Controlled Cracking. Nature, 485, pp.221-224.
10.1038/nature1100222575963
9
J.H., Kim and S. , Cho, 2014. Shape Design Sensitivity Analysis of Dynamic Crack Propagation Using Peridynamics and Parallel Computation. J. Comput. Struct. Eng. Inst. Korea, 27, pp.297-303.
10.7734/COSEIK.2014.27.4.297
페이지 상단으로 이동하기