Development of MLS Difference Method forMaterial Nonlinear Problem

Yoon Young-Cheol

doi:10.7734/COSEIK.2016.29.3.237

Preview

Journal of the Computational Structural Engineering Institute of Korea. 2016. 237-244
https://doi.org/10.7734/COSEIK.2016.29.3.237

Development of MLS Difference Method forMaterial Nonlinear Problem

MLS차분법을 이용한 재료비선형 문제 해석

Yoon Young-Cheol^†^*

Department of Civil Engineering, Myongji College, Seoul, 03656, Korea

•^*Corresponding author Yoon, Young-Cheol ycyoon@mjc.ac.kr

License (open-access, http://creativecommons.org/licenses/by-nc/3.0/):

This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

ABSTRACT

This paper presents a nonlinear Moving Least Squares(MLS) difference method for material nonlinearity problem. The MLS difference method, which employs strong formulation involving the fast derivative approximation, discretizes governing partial differential equation based on a node model. However, the conventional MLS difference method cannot explicitly handle constitutive equation since it solves solid mechanics problems by using the Navier’s equation that unifies unknowns into one variable, displacement. In this study, a double derivative approximation is devised to treat the constitutive equation of inelastic material in the framework of strong formulation; in fact, it manipulates the first order derivative approximation two times. The equilibrium equation described by the divergence of stress tensor is directly discretized and is linearized by the Newton method; as a result, an iterative procedure is developed to find convergent solution. Stresses and internal variables are calculated and updated by the return mapping algorithm. Effectiveness and stability of the iterative procedure is improved by using algorithmic tangent modulus. The consistency of the double derivative approximation was shown by the reproducing property test. Also, accuracy and stability of the procedure were verified by analyzing inelastic beam under incremental tensile loading.

본 연구는 재료비선형 문제를 다루기 위한 비선형 MLS 차분법의 정식화 과정을 제시한다. MLS 차분법은 절점모델을 기반으로 고속 미분근사식을 활용하여 지배 미분방정식을 직접 이산화 하는데, 변수를 변위로 일원화한 Navier 방정식을 사용하여 탄성재료 문제를 다룬 기존의 MLS 차분법은 재료의 구성방정식을 별도로 고려할 수 없다. 본 연구에서는 비선형 재료의 구성방정식을 반영할 수 있는 강정식화를 위해 1차 미분근사를 반복 사용하는 겹미분근사를 고안했다. 응력의 발산으로 표현되는 평형방정식을 그대로 이산화하고 Newton 방법을 적용하여 반복계산을 통해 수렴해를 찾는 비선형 알고리즘을 제시했다. 응력 계산과 내부변수의 갱신은 return mapping 알고리즘을 활용하였고, 알고리즘 접선계수(algorithmic tangent modulus)의 적용을 통해 빠르고 안정적인 반복계산이 가능하도록 하였다. 재생성 시험을 통해 겹미분근사의 정당성을 검증했고, 비선형재료에 대한 인장문제의 해석을 통해 개발된 비선형 MLS 차분 알고리즘의 정확성과 안정성을 확인하였다.

Keywords

material nonlinearity, MLS difference method, strong formulation, double derivative approximation, Newton method

재료비선형, MLS 차분법, 강정식화, 겹미분근사, Newton 방법

MAIN

1. 서 론

MLS(moving least squares) 차분법은 강정식화(strong formulation) 기반의 무요소법(meshfree method)으로서 요소망이나 그리드의 구성없이 Taylor 전개에 근거한 이동 최소제곱법을 사용한다. MLS 차분법은 최근까지 탄성균열 문제(Yoon et al., 2007), 복합재료 문제(Yoon and Lee, 2009), 동적해석 문제(Yoon et al., 2012), 동적균열전파 문제(Yoon et al., 2014) 등 다양한 분야에 적용되어 왔다. MLS 차분법의 근사함수는 요소를 사용하지 않기 때문에 요소의 경계를 따라 적합조건의 제약이 있는 유한요소법보다 편리하고, 재생성(reproducing property) 또는 일관성 (consistency)이 근사함수의 구성과정에서 자동으로 만족 되기 때문에 수학적 정당성도 자연스럽게 확보된다(Yoon and Song, 2014a; 2014b). 강정식화를 통해 약형식(weak form) 또는 변분식에 대한 적분이 필요없기 때문에 계산 효율성도 높다. 한편, 약정식화(weak formulation)를 기반으로 무요소 근사함수를 이용하여 재료비선형 문제를 다룬 연구들은 다수 있다(Gu et al., 2007; Pozo et al., 2009; Dai et al., 2006). 그러나 그 연구들은 모두 비선형 유한요소 해석 절차에서 근사함수만 무요소 근사함수로 대체하여 사용하는 방법으로 개발되었기 때문에 본 연구와 같이 강정식화를 기반으로 재료비선형성을 고려하는 데에는 활용이 어렵다.

본 연구에서는 재료비선형 문제를 다루기 위한 MLS 차분 알고리즘을 제시한다. 지금까지 MLS 차분법은 편의상 지배 미분방정식의 변수를 모두 변위로 바꾸어 Navier 방정식으로 불리는 미분방정식을 해석해 왔다. 그러나 Navier 방정식을 이산화하면 구성방정식이 외재적으로 나타나지 않기 때문에 기존에 개발된 재료비선형 모델을 적용할 수 없다. 따라서 본 연구에서는 응력텐서의 발산(divergence)으로 표현되는 지배 미분방정식을 그대로 이산화하면서 기존의 구성방정식을 활용 할 수 있는 해석법을 제시한다.

정적문제에서 체적력이 없는 경우 식 (1)의 Navier 방정식 에서는 $\nabla^{2} u$ 또는 ∇(∇·u) 와 같이 변위(u)의 2차 미분을 포함하는 고차항이 나타난다. 이 때, λ와 μ는 Lamé 상수 이다. 그러나 식 (2)와 같이 응력텐서(σ)의 발산으로 표현된 지배방정식에는 응력을 정의할 때 1차 미분항이 포함되고 발산을 계산할 때 다시 한번 1차 미분이 나타나지만 2차 미분항은 나타나지 않는다.

μ \nabla^{2} u + (λ + μ) \nabla (\nabla \cdot u) = 0

(1)

\nabla \cdot σ = 0

(2)

즉, 비선형 재료의 구성방정식을 고려하기 위해 2차 미분 대신 1차 미분을 두 번 계산하게 된다. 이와 같이 조금 어색한 정식화 절차로 인해 그 동안 강정식화를 사용하는 무요소법을 재료비선형 문제에 적용하려는 연구가 제대로 이루어지지 못했다.

본 연구의 주요 관심사는 이와 같은 미분근사 방법이 제대로 작동할 것인지의 여부와 함께 강정식화에 기반한 비선형 알고리즘이 제대로 작동하는지의 여부이다. 또한, Newton 방법 적용을 위한 잔차식의 구성, 기존 재료모델의 적용, 응력 갱신(update) 알고리즘도 함께 고려되어야 할 것이다. 재료비 선형성과 기하비선형성을 함께 고려하는 대변형 문제의 경우 유한요소법은 과도한 변형으로 인해 요소망이 심하게 찌그러 지면 참조(reference) 좌표계로의 매핑(mapping)이 어려워 지고 수치적분을 위한 적분점 배치가 곤란해지는 반면, 절점 계산만 수행하는 MLS 차분법은 이런 어려움이 없다. 그러나 지금까지 대변형 문제에 강정식화를 사용하는 무요소법이 적용된 사례는 없다.

본 논문에서 제시하는 MLS 차분법은 절점모델을 기반으로 이동최소제곱법을 통해 구성한 1차 미분근사식을 반복 사용 하여 응력텐서의 발산 형태로 주어진 지배 미분방정식을 그대로 이산화한다. 그 다음 비선형 재료에 대한 구성방정식을 적용 하고 반복계산을 통해 수렴해를 찾는다. 따라서 요소망 없이 전 해석과정을 진행하는 무요소법 본래의 모델링 장점이 극대화되고 수치적분이 없어 계산 효율성이 우수하다.

2. 1차원 재료비선형 문제의 정식화

2.1 1차원 지배방정식과 계방정식

미소변형을 가정한 1차원 고체역학 경계값 문제에서 체적 력을 무시할 때, 지배방정식과 경계조건은 다음과 같다.

\frac{\partial σ}{\partial x} = 0 in Ω

(3)

σ = \overline{t} in Γ_{t}

(4)

u = \overline{u} on Γ_{u}

(5)

여기서, $\overline{t}$ 는 표면력, $\overline{u}$ 는 변위경계값, ∂σ/∂x는 응력의 발산을 의미하고, Ω는 내부영역, $Γ_{t}$ 는 자연경계, $Γ_{u}$ 는 필수 경계를 가리킨다.

MLS 차분법의 미분근사식 유도를 위해 y를 기준으로 한 Taylor 급수로부터 m차 Taylor 다항식을 다음과 같이 정의 한다.

u_{L}^{m} (x, y) : = u (y) + (\frac{x - y}{ρ}) \frac{ρ}{1!} u^{(1)} (y) \dots

(6)

$+ {(\frac{x - y}{ρ})}^{m} \frac{ρ^{m}}{m!} u^{(m)} (y) = p_{m}^{T} (\frac{x - y}{ρ}) a (y)$

여기서, $p_{m}^{T} (\frac{x - y}{ρ})$ 는 m차 다항식 기저벡터이고, a(y)의 성분은 u(x)의 y 위치에서의 미분값들로 이루어진다. ρ는 이동최소제곱법 적용시 영향영역인데, 실제로 가중함수의 반경을 의미한다. 미분계수 u(y), u⁽¹⁾(y), …는 이동최소제곱 잔차식을 최소화하여 구한다.

선형 문제에서는 탄성계수가 변하지 않기 때문에 체적력이 없는 경우 Navier 방정식은 아래와 같다.

\frac{\partial σ (x)}{\partial x} = \frac{\partial}{\partial x} (E \frac{\partial u (x)}{\partial x}) = E \frac{{\partial u}^{2} (x)}{{\partial x}^{2}} = 0 i n Ω

(7)

그러나 재료비선형 문제에서 사용하는 탄소성 접선계수는 고정된 값이 아니고 응력상태에 따라 변하는 값이기 때문에 x에 대한 미분의 영향을 받아 식 (7)을 사용할 수 없다. 대신, 식 (3)의 미분방정식을 그대로 사용해야 구성방정식을 고려하고 기존의 비선형 재료모델도 사용할 수 있다.

따라서 본 연구에서는 탄소성 접선계수(E^E-P(x))가 외재적으로 나타나는 식 (8)을 정식화한다.

\frac{\partial}{\partial x} = (E^{E - P} (x) ε (x)) = \frac{\partial}{\partial x} (E^{E - P} \frac{\partial u (x)}{\partial x}) = 0

(8)

윗 식을 식 (7)과 비교하면 2차 미분항이 사라진 대신, 1차 미분항이 두 번 반복해서 나타나는 것을 알 수 있다. 즉, 1차 미분근사를 반복 사용하면 임의의 절점 x_J에 대해 다음과 같은 이산화가 가능하다.

\sum_{I}^{} Φ_{I}^{[(1)]} (x_{J}) σ_{I} = \sum_{I}^{} Φ_{I}^{[(1)]} (x_{J}) (E^{E - P} {(x}_{I}) ε_{I})

(9)

$= \sum_{I}^{} Φ_{I}^{[(1)]} (x_{J}) (E^{E - P} {(x}_{I}) \sum_{K}^{} Φ_{K}^{[(1)]} (x_{I}) u_{K})$

$= \sum_{I} \sum_{K}^{} Φ_{I}^{[(1)]} (x_{J}) (E^{E - P} {(x}_{I}) Φ_{K}^{[(1)]} (x_{I}) u_{K}$

여기서, $Φ_{I}^{[(1)]} (x_{J})$ 는 x_J에 대한 이웃절점인 절점 I에 대한 형상함수의 1차 미분근사이고, $Φ_{I}^{[(1)]} (x_{J})$ 는 그 반대이다. 흥미로운 점은 구성방정식 $E^{E - P} {(x}_{I}) ε_{I}$ 와 적합방정식 $ε_{I} = \sum_{K}^{} Φ_{K}^{[(1)]} (x_{I}) u_{K}$ 이 적용되면서 I에 대한 합(summation) 외에 K에 대한 합도 추가적으로 나타나는 것이다. 본 연구에서는 이와 같이 2차 미분근사를 위해 1차 미분근사를 두 번 반복하는 것을 겹미분근사(double derivative approximation)로 명명한다.

2차 미분근사 없이 탄소성 접선계수를 고려하기 위해 겹미분근사를 사용하여 내부영역의 절점들에 대해 식 (9)의 평형방정식을 조립(assemble)하면 다음과 같은 계방정식 (system of equations)을 얻는다.

A_{J}^{N} = 1 (\sum_{I} Φ_{I}^{[(1)]} (x_{J}) σ_{I}) = (Φ^{[1]}) D Ψ^{[1]} Δu

(10)

여기서, $A_{J = 1}^{N}$ 은 J=1,···,N에 대한 조립을 의미한다. 식 (10)에 나타난 각 항들은 다음과 같이 정의된다.

D = Diag (E^{E - P} (x_{1}), E^{E - P} (x_{2}), \dots, E^{E - P} (x_{N}))

(11)

Ψ^{[1]} = (\begin{matrix} Φ_{1}^{[(1)]} (x_{1}) & \dots & \begin{matrix} Φ_{N}^{[(1)]} (x_{1}) \end{matrix} \\ . . . & . . . & . . . \\ \begin{matrix} Φ_{1}^{[(1)]} (x_{N}) \end{matrix} & \dots & \begin{matrix} Φ_{N}^{[(1)]} (x_{N}) \end{matrix} \end{matrix})

(12)

Δu = {({Δu}_{1}, \dots, {Δu}_{N})}^{T}

(13)

위의 식들은 모두 전체(global) 좌표계에서 만들어졌다는 것을 주목할 필요가 있다. 국부(local) 좌표계를 사용하지 않기 때문에 프로그래밍 과정이 단순하다. Ψ_[1]는 x_I에 대한 이웃절점들의 1차 미분근사값을 절점순서에 따라 행별로 내부절점들에 대해 조립한 것에 불과하다. 즉, Ψ_[1]의 열은 $Φ_{I}^{[(1)]} (x_{J})$ 값들을 포함하고, 행은 $Φ_{K}^{[(1)]} (x_{I})$ 를 포함한다. 한편, 식 (14)에서 보듯이, 식 (9)의 평형방정식은 식 (10)의 계방정식의 행벡터와 동일한 의미를 갖는다.

{\sum_{I} Φ}_{I}^{[(1)]} (x_{J}) σ_{I} = ({Row}^{K - t h} Ψ^{[1]}) {DΨ}^{[1]} Δu

(14)

결과적으로 2차 미분근사를 사용하는 기존의 강정식화에 비해 겹미분근사를 사용한 새로운 강정식화는 지배방정식의 이산화 과정을 복잡하게 만들지도 않고, 계산량을 증가 시키지도 않는다. 1차원 문제의 정식화 과정에서 Γ_t 와 Γ_u 에 대한 경계조건은 기존과 동일한 방법으로 이산화되므로 다음 절에서 설명하는 반복계산 절차 외에 재료비선형성 고려로 인한 추가적인 계산은 발생하지 않는다.

2.2 반복계산을 위한 지배방정식의 선형화

재료비선형 문제를 풀기 위해서는 반복계산을 통해 수렴해를 찾는 알고리즘이 필요하다. 식 (15)는 지배 미분방정식을 Newton 방법으로 선형화시킨 잔차식(residual equation)으로, i+1번째 단계(step)에서 (k+1)번째 반복계산(iteration)을 수행하기 위한 잔차식이다.

\hat{R}_{i + 1}^{(k + 1)} (x) = \hat{R}_{i + 1}^{(k)} (x) + \frac{\partial}{{\partial u}_{i + 1}^{(k)}} (\frac{\partial σ_{i + 1}^{(k)}}{\partial x}) Δ u_{i + 1}^{(k + 1)}

(15)

$= \hat{R}_{i + 1}^{(k)} (x) + \frac{\partial}{\partial x} (\frac{\partial σ_{i + 1}^{(k)}}{\partial ε_{i + 1}^{(k)}} \frac{\partial ε_{i + 1}^{(k)}}{\partial u_{i + 1}^{(k)}}) Δ u_{i + 1}^{(k + 1)}$

여기서, $Δ u_{i + 1}^{(k + 1)}$ 은 현재 반복계산에서 구하고자 하는 변위 증분이고, $\hat{R}_{i + 1}^{(k)} (x) = \partial σ_{i + 1}^{(k)} / \partial x$ 은 (k)번째 반복계산시 응력의 발산에 대한 잔차이다. $\partial σ_{i + 1}^{(k)} / \partial ε_{i + 1}^{(k)}$ 은 탄소성 접선계수를 가리키며, 다차원 문제의 경우 계산의 효율성과 정확성 확보를 위해 알고리즘 접선계수(algorithmic tangent modulus)를 사용한다. 식 (15)에서 u에 대한 미분과 에 대한 미분의 위치를 교환한 것을 주목할 필요가 있다. $σ_{i + 1}^{(k)} (x)$ 은 return mapping 알고리즘(Simo and Taylor, 1985)으로 계산되며, 역학적(kinematic) 변수들은 backward Euler 방법으로 갱신된다(Simo and Hughes, 1998).

$\hat{R}_{i + 1}^{(k + 1)} (x) = 0$ 을 가정하고, 식 (10)의 계방정식을 고려 하면, x_J에서 계산한 식 (15)는 다음과 같다.

\hat{R}_{i + 1}^{(k)} (x_{J}) = - ({Row}^{J - t h} Ψ^{[1]}) D_{i + 1}^{(k)} Ψ^{[1]} {Δu}_{i + 1}^{(k + 1)}

(16)

여기서, Row_J-th(·)는 ·의 J번째 행벡터를 가리킨다. 위의 잔차식을 내부영역 절점들에 대해 조립한 후 변위증분에 대해 풀면 다음과 같은 전체 계방정식을 얻는다.

{Δu}_{i + 1}^{(k + 1)} = - {(Ψ^{[1]} D_{i + 1}^{(k)} Ψ^{[1]})}^{- 1} R_{i + 1}^{(k)}

(17)

여기서, $R_{i + 1}^{(k + 1)} = {(\hat{R}_{i + 1}^{(k)} (x_{1}), \dots, \hat{R}_{i + 1}^{(k)} (x_{N}))}^{T}$ 은 내부영역 모든 절점에 대해 $\hat{R}_{i + 1}^{(k)} (x_{J})$ 값을 조립한 전체 벡터이다.

결과적으로 1차 겹미분근사를 이용하여 탄소성 접선계수를 외재적으로 끌어내고 재료비선형성을 갖는 2차 미분방정식을 강정식화 하였다. 반복계산의 정확성과 안정성이 확보되면 고차미분 계산에 대한 부담이 없기 때문에 재료비선형 문제의 효율적 해석이 가능하며 다차원 문제에도 그대로 확장 가능하다.

3. 다차원 재료비선형 문제의 정식화

3.1 고속 미분근사

Taylor 전개를 기반으로 이동최소제곱법을 사용하여 형상 함수를 도출하는 MLS 차분법의 미분근사식 유도과정은 이미 여러 편의 논문들에서 소개된 바 있다(Yoon and Lee, 2004; 2009; Yoon et al.,, 2007; 2012; 2014). 본 절에서는 다중지수 표기법(multi-index notation)을 사용하여 미분 근사식을 간단히 소개한다. x=(x₁,…,x_n)는 n차원 좌표를 갖는 벡터이고, x의 α-차 지수와 편미분은 각각 $x^{α} = x_{1}^{α_{1}}$ $x_{2}^{α_{2}}, \dots, x_{n}^{α_{n}}$ 와 $D_{x}^{α} = \partial_{x_{1}}^{α_{1}} \partial_{x_{2}}^{α_{2}} \dots \partial_{x_{n}}^{α_{n}}$ 로 표현된다. m차 이상의 고차항을 무시하면, u(x,y)에 대한 m차 Taylor 다항식은 식 (18)과 같다.

u^{m} (x,y) = (\frac{{(x - y)}^{α_{1}}}{α_{1}!}, \dots, \frac{{(x - y)}^{α_{K}}}{α_{K}!}) (\begin{matrix} D_{x}^{α_{1}} u (y) \\ . . . \\ \begin{matrix} D_{x}^{α_{K}} u (y) \end{matrix} \end{matrix})

(18)

$= p_{m}^{T} (x, y) c (y)$

여기서, $K = \frac{(m + n)!}{m! n!}$ , α₁=(0,…,0), α_K(0,…,)이고, $= p_{m}^{T} (x, x_{0})$ 는 m차 다항식 기저벡터, c(y)는 $D_{x}^{α} (u (y))$ 의 각 성분들을 순서대로 정렬한 미지벡터이다. 이동최소제곱법 절차에 따라 절점별 가중치를 고려한 절점해의 잔차식을 구성하고 그것을 최소화하면 다음과 같이 c(y)를 형상함수의 근사미분과 절점해의 곱으로 얻는다.

c (y) = (\begin{matrix} D_{x}^{α_{1}} u (y) \\ . . . \\ \begin{matrix} D_{x}^{α_{K}} u (y) \end{matrix} \end{matrix}) ≅ (\begin{matrix} \begin{matrix} ∅_{1}^{[α_{1}]} u (y) \end{matrix} & \dots & \begin{matrix} \begin{matrix} ∅_{N}^{[α_{1}]} u (y) \end{matrix} \end{matrix} \\ . . . & . . . & . . . \\ \begin{matrix} \begin{matrix} ∅_{1}^{[α_{K}]} u (y) \end{matrix} \end{matrix} & \dots & \begin{matrix} \begin{matrix} ∅_{N}^{[α_{K}]} u (y) \end{matrix} \end{matrix} \end{matrix}) \cdot (\begin{matrix} u_{1} \\ . . . \\ u_{N} \end{matrix})

(19)

여기서, u_I(=u(x_I))는 x_I에서의 절점해이고, $∅_{I}^{[α]} (x)$ 는 절점 I에 대한 형상함수의 α차 미분근사를 의미한다. 식 (18)에 식 (19)를 대입하면서 y를 x로 치환하면 최종적인 Taylor 다항식이 얻어진다. MLS 차분법은 c(x)의 계산시 α_K차 까지의 모든 미분근사식을 한꺼번에 계산하기 때문에 계산 속도가 빠르고, 기저벡터의 차수 조정을 통해 미분근사의 정밀도 조정도 쉽다.

3.2 Newton 방법을 이용한 비선형 해석

MLS 차분법은 재료비선형 문제의 정식화시 탄소성 접선 계수의 적용을 위해 응력텐서의 발산 형태로 주어지는 평형 방정식을 그대로 이산화하며 이 과정에서 1차 겹미분근사를 사용한다. 체적력이 없을 때, 평형방정식은 Newton 방법을 적용하여 변위에 대한 미분과 그 증분을 이용하여 식 (20)과 같이 쓸 수 있다.

\frac{\partial}{\partial u} (\nabla \cdot σ) Δu = \frac{\partial}{\partial u} (\frac{\partial σ}{\partial x}) Δu = \frac{\partial}{\partial x} (\frac{\partial σ}{\partial ε} \frac{\partial ε}{\partial u}) Δu

(20)

여기서, 변위(u)에 대한 미분과 위치(x)에 대한 미분을 서로 교환한 것을 주목할 필요가 있다. 이것은 일반적인 수학식의 전개에서 항상 적용되지는 않지만 MLS 차분법의 미분 근사식을 이용하여 이산화할 때 미분의 순서가 계산결과에 영향을 주지 않는다는 점을 이용한 것이다. 또한, $\partial / \partial x$ 와 $\partial ε / \partial u$ 는 이산화 과정에서 유한요소법의 ‘B’ 행렬과 유사한 형태로 나타나기 때문에 식 (20)이 강정식화 임에도 불구하고 대칭성을 가질 것으로 예측된다.

식 (20)의 평형방정식을 반복계산을 위해 선형화하면 다음과 같다.

\hat{R}_{i + 1}^{(k + 1)} (x) = \hat{R}_{i + 1}^{(k)} (x) + \frac{{\partial σ}_{i + 1}^{(k)}}{{\partial u}_{i + 1}^{(k)}} Δ \hat{u}_{i + 1}^{(k + 1)}

(21)

$= \hat{R}_{i + 1}^{(k)} (x) + \frac{\partial}{\partial x} (\frac{{\partial σ}_{i + 1}^{(k)}}{{\partial ε}_{i + 1}^{(k)}} \frac{{\partial ε}_{i + 1}^{(k)}}{{\partial u}_{i + 1}^{(k)}}) Δ \hat{u}_{i + 1}^{(k + 1)}$

윗 식에서 ${\partial σ}_{i + 1}^{(k)} / {\partial ε}_{i + 1}^{(k)}$ 은 탄소성 접선계수를 의미한다. 또한, $\hat{R}_{i + 1}^{(k)} (x)$ 의 계산시 필요한 $σ_{i + 1}^{(k)}$ 은 탄성예측(elastic prediction)과 소성보정(plastic correction)을 이용한 return mapping 알고리즘으로 계산된다(Simo and Taylor, 1985). 다양한 역학적 변수들은 backward Euler방법으로 갱신한다 (Simo and Hughes, 1998).

이제 $\hat{R}_{i + 1}^{(k)} (x)$ =0을 가정하고 식 (21)을 MLS 차분법의 미분근사식으로 이산화하면 다음 식을 얻는다.

\hat{R}_{i + 1}^{(k)} (x_{J}) = - \frac{\partial}{\partial x} (\frac{{\partial σ}_{i + 1}^{(k)}}{{\partial ε}_{i + 1}^{(k)}} \frac{{\partial ε}_{i + 1}^{(k)}}{{\partial u}_{i + 1}^{(k)}}) Δ \hat{u}_{i + 1}^{(k + 1)}

(22)

$= - {Row}^{J - t h} {(Ψ^{[1]})}^{T} D_{i + 1}^{[k]} Ψ^{[1]} Δ u_{i + 1}^{(k + 1)}$

여기서, $\hat{R}_{i + 1}^{(k)} (x_{J})$ 는 직전 반복계산에서 return mapping 알고리즘으로 얻은 응력의 발산을 계산한 벡터이고, ${Δu}_{i + 1}^{(k + 1)}$ 는 모든 절점의 변위증분을 포함한 전체 벡터이다.

내부절점에 대해 식 (22)의 잔차식을 조립하고 변위증분에 대해 정리하면 다음과 같은 계방정식을 얻는다.

{Δu}_{i + 1}^{(k + 1)} = - {[{(Ψ^{[1]})}^{T} D_{i + 1}^{(k)} Ψ^{[1]}]}^{- 1} R_{i + 1}^{{k}}

(23)

여기서, $R_{i + 1}^{{k}}$ 는 내부영역 모든 절점에 대해 ${\hat{R}}_{i + 1}^{{k}} (x_{J})$ 값들을 조립한 전체 벡터이고, 따라서 ${\hat{R}}_{i + 1}^{{k}} (x_{J})$ 는 $R_{i + 1}^{{k}}$ 의 J번째 성분 벡터가 된다. 전체 계방정식 유도과정에서 1차원 문제의 경우 Ψ_[1]의 전치행렬이 필요없으나 다차원 문제는 $Ψ^{{[1]}^{T}}$ 가 나타나는 것을 주의할 필요가 있다. $D_{i + 1}^{(k)}$ 는 $D_{i + 1}^{(k)} (x_{1}), \dots, D_{i + 1}^{(k)} (x_{N})$ 를 대각선 성분으로 갖는 식 (24)와 같은 대각행렬 이다.

D_{i + 1}^{(k)} = [\begin{matrix} {\hat{D}}_{i + 1}^{(k)} (x_{1}) \\ 0 & . . . & 0 \\ {\hat{D}}_{i + 1}^{(k)} (x_{N}) \end{matrix}]

(24)

이 때, $D_{i + 1}^{(k)}$ 는 미분근사식에서 다중지수로 표기된 미분기호 $D_{x}^{α}$ 와 구분되어야 한다. 2차원 문제를 가정하면, x_I에 대한 재료행렬식 ${\hat{D}}_{i + 1}^{(k)} (x_{I})$ 는 다음과 같다.

(25)

다음으로 $Ψ^{[1]}$ 와 그 성분을 이루고 있는 $Ψ_{IJ}^{[1]}$ 은 식 (26)과 식 (27)에 제시했다.

Ψ_{IJ}^{[1]} = [\begin{matrix} Φ_{11}^{[1]} & \dots & \begin{matrix} Φ_{1 N}^{[1]} \end{matrix} \\ . . . & . . . & . . . \\ \begin{matrix} Φ_{\begin{matrix} \begin{matrix} N \end{matrix} \end{matrix} 1}^{[1]} \end{matrix} & \dots & \begin{matrix} Φ_{\begin{matrix} \begin{matrix} N \end{matrix} \end{matrix} \begin{matrix} \begin{matrix} N \end{matrix} \end{matrix}}^{[1]} \end{matrix} \end{matrix}]

(26)

{\begin{matrix} Φ \end{matrix}}_{IJ}^{[1]} = [\begin{matrix} Φ_{J}^{[(1, 0)]} (x_{I}) & \begin{matrix} 0 \end{matrix} \\ 0 & \begin{matrix} Φ_{J}^{[(0, 1)]} (x_{I}) \end{matrix} \\ \begin{matrix} \begin{matrix} Φ_{J}^{[(0, 1)]} (x_{I}) \end{matrix} \end{matrix} & \begin{matrix} \begin{matrix} Φ_{J}^{[(1, 0)]} (x_{I}) \end{matrix} \end{matrix} \end{matrix}]

(27)

여기서, ${\begin{matrix} Φ \end{matrix}}_{IJ}^{[1]}$ 의 I, J 아래첨자가 절점 I(또는 J)의 이웃절점에 포함된 절점 J(또는 x_J)의 형상함수와 연관된다는 점을 주목할 필요가 있다. 식 (26)의 전치행렬은 다음과 같다.

{(Ψ^{[1]})}^{T} = [\begin{matrix} {Φ_{11}^{[1]}}^{T} & \dots & \begin{matrix} {Φ_{1 N}^{[1]}}^{T} \end{matrix} \\ . . . & . . . & . . . \\ \begin{matrix} {Φ_{\begin{matrix} \begin{matrix} N \end{matrix} \end{matrix} 1}^{[1]}}^{T} \end{matrix} & \dots & \begin{matrix} {Φ_{\begin{matrix} \begin{matrix} N \end{matrix} \end{matrix} \begin{matrix} \begin{matrix} N \end{matrix} \end{matrix}}^{[1]}}^{T} \end{matrix} \end{matrix}]

(28)

따라서 식 (22)의 ${Row}^{J - t h {(Ψ^{[1]})}^{T}}$ 는 식 (26)으로부터 다음과 같이 표현된다.

{Row}^{J - t h {(Ψ^{[1]})}^{T}} = [{(Φ_{J1}^{[1]})}^{T}, \dots, {(Φ_{JN}^{[1]})}^{T}]

(29)

유한요소법은 요소별로 적분에 따른 국부 행렬식을 구성한 후 조립하지만, MLS 차분법은 적분 과정이 필요없고 실제로 국부 행렬식에 대한 계산이 없기때문에 처음부터 전체 계방 정식을 조립할 수 있는 장점이 있다. 이것은 전체 계방정식 식 (23)을 처음부터 만들어 나갈 수 있다는 의미이다. 만약, 절점을 유한요소법의 요소처럼 취급한다면, 요소의 평형 방정식에 해당하는 MLS 차분법의 절점에서의 평형방정식은 식 (21)로부터 다음과 같이 쓸 수 있다.

\frac{\partial}{\partial x} [\frac{\partial σ_{i + 1}^{(k)}}{\partial ε_{i + 1}^{(k)}} \frac{\partial ε_{i + 1}^{(k)}}{\partial u_{i + 1}^{(k)}}] Δ {\hat{u}}_{i + 1}^{(k + 1)}

(30)

윗 식을 보면 식 (20)의 $\partial / \partial x$ 와 $\partialε / \partial u$ 가 이산화 과정에서 각각 ${(Φ_{IJ}^{[1]})}^{T}$ 와 $Φ_{IJ}^{[1]}$ 로 나타난 것을 알 수 있다.

자연경계조건과 필수경계조건은 기존과 같은 방식으로 이산화 할 수 있으므로 식 (23)의 계방정식을 조립하는 과정 에서 자연경계와 필수경계에 배치된 절점들에 대한 잔차식을 함께 포함시키면 전체 계방정식이 완성된다. 이렇게 전체 계방정식이 구성되면 본 절에서 제시된 Newton 방법의 절차에 따라 반복계산을 통해 수렴해를 찾을 수 있다. 좀 더 자세한 과정은 Simo와 Hughes(1998)를 참고할 수 있다.

4. 수치예제

4.1 1차 겹미분근사의 정확성 검증

본 절에서는 1차 미분을 두 번 사용하여 2차 미분을 근사 하는 1차 겹미분근사의 정확성을 검증한다. 이론해를 알고 있는 간단한 1차원 탄성봉 문제를 고려한다. 기지의 변위를 이용하여 기존의 2차 미분근사와 1차 겹미분근사를 사용하여 응력의 발산 값을 재생하고 오차를 측정하여 미분근사의 정확성을 검증한다. 이와 같은 과정을 재생성 시험으로 부르기도 한다 (Lee and Yoon, 2004; Yoon and Song, 2014b).

Fig. 1은 1차원 봉 문제에 대한 재생성 시험결과를 보여 준다. 기존의 2차 미분근사와 비교할 때 1차 겹미분근사의 수렴성은 동일한 기울기를 나타내고 절대적인 오차의 크기는 오히려 약간 작게 나타났다. 이것은 1차 겹미분근사가 2차 미분근사와 동등하거나 조금 더 나은 성능을 보인다는 것을 의미한다. 따라서 1차 겹미분근사를 재료비선형 문제의 해석을 위한 강정식화에 적용하는 정당성을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2016-029-03/04TK062016290304/images/10.7734.29.3.237.F100.png

Figure 1

Reproducing property of 1st order double derivative approximation

4.2 1차원 비탄성 봉의 인장문제

본 절에서는 비선형 MLS 차분 알고리즘을 이용하여 1차원 비탄성 봉의 인장문제를 해석하였다. 해석에 사용한 물성치는 E=1×10⁶, σ_y=1×10⁴, H=5×10³(경화계수), K=5×10⁴(소성계수)이다. 61개의 절점과 2차 다항식을 사용했다.

Fig. 2는 변위와 반력의 관계인데 선형강화 재료에 대한 응답거동을 보여주며, Fig. 3은 해석모델의 중앙에 위치한 절점에서 계산한 응력과 변형률 관계인데 재료의 항복 이후 선형강화 거동을 정확하게 모사하고 있다.

해석결과, 재료의 항복이후 소성영역에서 반복계산의 횟수는 모두 5회 미만으로 나타나 비선형 MLS 차분 알고리즘의 수렴성이 우수함을 확인하였다. 반복계산을 통해 수렴해를 찾는 해석 알고리즘의 안정성을 추가적으로 조사하기 위해 매 해석단계에서 계산된 변위값에 대한 L² 상대오차를 해석 단계별로 Fig. 4에 도시했다. 모든 단계에서 오차의 절대값은 1.0×10^-10이하의 크기로 나타났으며, 재료가 소성영역에 들어가는 200단계 이후에도 상당히 안정적인 오차 거동을 보이는 것으로 보아 비선형 MLS 차분법이 효과적으로 작동하는 것을 알 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2016-029-03/04TK062016290304/images/10.7734.29.3.237.F101.png

Figure 2

Displacement and reaction curve(1-D rod)

https://cdn.apub.kr/journalsite/sites/jcoseik/2016-029-03/04TK062016290304/images/10.7734.29.3.237.F102.png

Figure 3

Stress and strain response at a center node(1-D rod)

https://cdn.apub.kr/journalsite/sites/jcoseik/2016-029-03/04TK062016290304/images/10.7734.29.3.237.F103.png

Figure 4

Relative L² error norm according to loading steps(1-D rod)

4.3 2차원 비탄성 보의 인장문제

본 절에서는 앞 절에서 다루었던 1차원 비탄성 봉의 인장문제를 Fig. 5과 Fig. 6에서 보듯이 2차원 모델로 모형화하여 해석하였다. 해석에 사용한 물성치는 1차원 봉 문제와 동일하고, 153(=17×9)개의 절점과 2차 다항식을 사용했다. 1차원 문제와의 비교를 위해 포아송비를 0으로 적용했으며, 평면응력 조건을 고려하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2016-029-03/04TK062016290304/images/10.7734.29.3.237.F104.png

Figure 5

2-D inelastic beam problem under uniform tension load

https://cdn.apub.kr/journalsite/sites/jcoseik/2016-029-03/04TK062016290304/images/10.7734.29.3.237.F105.png

Figure 6

Node distribution of 2-D inelastic beam problem

Fig. 7에는 2차원 모델로 해석한 보의 변위와 반력의 관계를 1차원 모델 해석결과와 함께 도시했다. 2차원 해석 결과로 얻은 탄성구간과 변형률경화 구간의 응답이 1차원 경우의 응답과 정확하게 일치하는 것을 확인할 수 있다. Fig. 8은 2차원 모델의 중앙에 위치한 절점에서 계산한 응력과 변형률의 관계인데, 역시 2차원 해석결과가 1차원 해석결과와 잘 일치하는 것을 볼 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2016-029-03/04TK062016290304/images/10.7734.29.3.237.F106.png

Figure 7

Displacement and reaction curve(2-D beam)

https://cdn.apub.kr/journalsite/sites/jcoseik/2016-029-03/04TK062016290304/images/10.7734.29.3.237.F107.png

Figure 8

Stress and strain response at a center node (2-D beam)

2차원 보 문제의 경우 소성영역에서 반복계산의 횟수는 5~6회 정도로 나타났다. 알고리즘 접선계수를 사용하여 계산 효율성을 확보하였고, 탄성에서 소성으로 넘어갈 때에도 갑작스런 반복계산 횟수의 증가없이 안정적인 수렴성을 보였다. 모든 해석단계에서 반복계산이 10회를 넘지 않았는데 이것은 1차 겹미분근사가 재료비선형 문제의 강정식화가 필요로 하는 미분계산의 정확도를 제공하는 것을 의미한다. 결국, 2차 미분이 포함된 지배방정식을 약정식화 없이도 1차 미분근사 만으로 정확하게 해석할 수 있음을 보여준다.

5. 결 론

본 연구에서는 수치적분이나 필수경계조건 처리를 위한 구속조건식의 도입없이 무요소법의 장점을 최적화할 수 있는 강정식화에 기반한 MLS 차분법을 이용하여 재료비선형 문제를 해석할 수 있는 알고리즘을 제시하였다. MLS 차분법은 요소나 그리드에 의존적이지 않고 기존 무요소법보다 빠르게 미분 근사를 하기 때문에 강정식화에 유리하다. 기존의 MLS 차분법은 지배방정식을 변위만의 식인 Navier 방정식으로 변환하여 사용했으나, 본 연구에서는 비선형 재료의 구성방정식을 활용 하기 위해 응력텐서의 발산으로 표현되는 평형방정식을 그대로 이산화하였다. 따라서 응력의 발산값을 계산하는 과정과 변위의 미분항이 포함된 응력을 계산하는 과정에서 1차 미분계산이 두 번 필요하다. 본 연구는 절점만 사용하는 MLS 차분법의 유연성을 활용하여 약정식화의 도입없이 1차 미분근사를 반복하는 겹미분근사를 고안하여 강정식화 기반의 계방정식을 유도하고, Newton 방법에 근거하여 반복계산이 가능한 비선형 알고리즘을 제시했다. 또한, 응력값의 계산과 역학적 변수들의 갱신은 return mapping 알고리즘을 이용하였고, 알고리즘 접선계수를 적용하여 반복계산의 안정성과 효율성을 높였다.

겹미분근사의 수학적 검증을 위해 재생성 시험을 수행하고 1차 겹미분근사가 기존의 2차 미분근사와 비교하여 동등 하거나 더 나은 정확성을 갖는 것을 확인했다. 개발된 비선형 알고리즘의 정확성과 안정성을 확인하기 위해 1차원 비탄성 봉의 인장문제를 해석하여 선형강화 재료에 대한 정확한 응답을 얻을 수 있음을 보였다. 2차원 보의 인장 문제에서도 재료의 비선형성을 정확히 반영하는 변위, 변형률, 응력을 얻었고, 1차원 문제의 해석결과와도 잘 일치하는 것을 확인 했다. 수렴해를 찾기 위한 반복계산은 1차원 문제에서 5회 미만, 2차원 문제에서는 10회 미만으로 나타나 겹미분근사를 활용한 비선형 알고리즘이 알고리즘 접선계수와 함께 잘 작동 하는 것을 보였다. 결과적으로 개발된 비선형 MLS 차분법이 재료비선형 문제를 정확하게 해석할 수 있음을 확인하였으며, 향후 재료비선형과 기하비선형을 모두 고려한 대변형 문제의 해석까지 확장될 것으로 기대된다.

Acknowledgements

이 논문은 2014년도 정부(미래창조과학부)의 재원으로 한국연구재단의 지원을 받아 수행된 기초연구사업임(No. 2014R1A1A1002000).

References

K.Y., Dai, G.R., Liu, X., Han and Y., Li, 2006. Inelastic Analysis of 2D Solids using a Weak-form RPIM based on Deformation Theory. Comput. Methods Appl. Mech. & Eng., 195, pp.4179-4193.

10.1016/j.cma.2005.07.021

Y.T., Gu, Q.X., Wang, K.Y., Lam and K.Y., Dai, 2007. A Pseudo-elastic Local Meshless Method for Analysis of Material Nonlinear Problems in Solids. Eng. Anal. Bound. Elem., 31, pp.771-782.

10.1016/j.enganabound.2006.12.008

P.L., Pozo, F., Perazzo and A., Angulo, 2009. A Meshless FPM Model for Solving Nonlinear Material Problems with Proportional Loading based on Deformation Theory. Adv. Eng. Softw., 40, pp.1148-1154.

10.1016/j.advengsoft.2009.05.004

S.H., Lee and Y.C., Yoon, 2004. Meshfree Point Collocation Method for Elasticity and Crack Problems. Int. J. Numer. Methods Eng., 61, pp.22-48.

10.1002/nme.1053

J.C., Simo and R.L., Taylor, 1985. Consistent Tangent Operators for Rate-independent Elastoplasticity. Comput. Methods Appl. Mech. & Eng., 48, pp.101-118.

10.1016/0045-7825(85)90070-2

Simo J.C. Hughes T.J.R. 1998Computational inelasticitySpringer-VerlagNew York

Y.C., Yoon, D.J., Kim and S.H., Lee, 2007. A Gridless Finite Difference Method for Elastic Crack Analysis. J. Comput Struct. Eng., 20, pp.321-327.

Y.C., Yoon and S.H., Lee, 2009. Intrinsically Extended Moving Least Squares Finite Difference Method for Potential Problems with Interfacial Boundary. J. Comput Struct. Eng., 22, pp.411-420.

Y.C., Yoon, K.H., Kim and S.H., Lee, 2012. Dynamic Algorithm for Solid Problems using MLS Difference Method. J. Comput. Struct. Eng., 25, pp.139-148.

10.7734/COSEIK.2012.25.2.139

Y.C., Yoon, K.H., Kim and S.H., Lee, 2014. Analysis of Dynamic Crack Propagation using MLS Difference Method. J. Comput. Struct. Eng., 27, pp.17-26.

10.7734/COSEIK.2014.27.1.17

Y.C., Yoon and J.-H., Song, 2014a. Extended Particle Difference Method for Weak and Strong Discontinuity Problems: Part I. Derivation of the Extended Particle Derivative Approximation for the Representation of Weak and Strong Discontinuities. Comput. Mech., 53, pp.1087-1103.

10.1007/s00466-013-0950-8

Y.C., Yoon and J.-H., Song, 2014b. Extended Particle Difference Method for Weak and Strong Discontinuity Problems: Part II. Formulations and Applications for Various Interfacial Singularity Problems. Comput. Mech., 53, pp.1105-1128.

10.1007/s00466-013-0951-7

Journal of the Computational Structural Engineering Institute of Korea ISSN:1229-3059(Print) 2287-2302(Online) 한국전산구조공학회 논문집

Preview

Development of MLS Difference Method forMaterial Nonlinear Problem

ABSTRACT

MAIN

Figure 1

Reproducing property of 1st order double derivative approximation

Figure 2

Displacement and reaction curve(1-D rod)

Figure 3

Stress and strain response at a center node(1-D rod)

Figure 4

Relative L2 error norm according to loading steps(1-D rod)

Figure 5

2-D inelastic beam problem under uniform tension load

Figure 6

Node distribution of 2-D inelastic beam problem

Figure 7

Displacement and reaction curve(2-D beam)

Figure 8

Stress and strain response at a center node (2-D beam)

Acknowledgements

References

Relative L² error norm according to loading steps(1-D rod)