Development of Meshless Method Using Least-Squares Method with Geometric Conservation Law for Structural Dynamic Analysis

Sang Woo Lee; Jin Young Huh; Kyu Hong Kim

doi:10.7734/COSEIK.2023.36.1.67

Preview

Research Paper

Journal of the Computational Structural Engineering Institute of Korea. 28 February 2023. 67-74
https://doi.org/10.7734/COSEIK.2023.36.1.67

Development of Meshless Method Using Least-Squares Method with Geometric Conservation Law for Structural Dynamic Analysis

기하학적 보존을 만족하는 최소제곱법을 활용한 무격자 구조해석 기법 개발

Sang Woo Lee¹^*

Jin Young Huh²

Kyu Hong Kim³

이 상우¹^*

허 진영²

김 규홍³

¹Graduate Student, Aerospace Engineering, Seoul National University, Seoul, 08826, Korea

²Senior Researcher, Agency for Defense Development, Daejeon, 34186, Korea

³Professor, Aerospace Engineering, Seoul National University/Institution of Advanced Aerospace Technology, 08826, Korea

¹서울대학교 항공우주공학과 박사과정

²국방과학연구소 선임연구원

³서울대학교 항공우주공학과/항공우주신기술연구소 교수

^{*Corresponding Author}

ABSTRACT

A meshless technique using the geometric conservation least-squares method (GC-LSM) was devised to discretize the governing equation of linear elasticity. Although the finite-element method is widely used for structural analysis, a meshless method was developed because of its advantages in a moving grid system. This work is the preliminary phase for developing a fully meshless-based fluid-structure interaction solver. In this study, Cauchy’s momentum equation was discretized in strong form using GC-LSM for the structural domain, and the Newmark beta method was used for time integration. The solver was validated in 1D, 2D, and 3D benchmarking problems. Static and dynamic results were obtained. The results are more accurate than those of analytic solutions.

Keywords

computational structure dynamics

meshless method

geometric conservation least-squares method (GC-LSM)

dynamic analysis

strong form

본 논문은 코시 모멘텀(Cauchy’s Momentum) 방정식을 이산화하기 위해 GC-LSM(Geometric Conservative Least Squares Method)을 사용한 새로운 Meshless 방법을 제시한다. FEM(Finite Element Method) 방법이 구조해석에 널리 사용되고 있지만 무격자 기법은 격자를 이동해야 할 때 장점이 많기 때문에 개발되었다. 본 작업은 무격자 기반의 FSI(Fluid-Structure Interaction) 프로그램을 개발하기 위한 기틀을 다지는 단계이다. 본 논문에서는 Cauchy’s Momentum 방정식을 GC-LSM을 사용하여 강형식 형태로 이산화하였고, 시간 적분을 위해 New Mark Beta 방법을 사용하였다. 개발된 기법은 1D, 2D 및 3D 벤치마킹 문제에서 검증했으며, 정적 해석 및 동적 해석 결과가 해석해와 비교시 매우 정확한 결과를 보여준다.

키워드

전산구조해석

무격자 기법

동적해석

강형식

기히학적 보존 최소제곱법

MAIN

1. 서 론
2. 수치기법
2.1 지배방정식
2.2 구성방정식
2.3 무격자 기법
2.4 시간전진 기법
2.5 경계조건
2.4 수치해석 결과 및 고찰
3. 결 론

1. 서 론

최근 개발되는 고성능 전투기나 민항기의 경우 탄성을 가진 구조물이 큰 비율을 차지하며, 공탄성에 의한 영향이 두드러진다(Slotnick et al., 2014). 이러한 공탄성의 영향으로 날개가 심하게 진동하거나, 대변형을 일으켜 비행체의 파손이 발생할 수도 있다. 또한 비행체가 파손이 되지 않더라도 날개의 진동현상으로 공력 특성이 변해 비행 성능이 설계 목표보다 저하될 수 있고, 진동이 비행체 내부로 전해져 전자장비나 탑재체 등의 악영향을 끼칠 수 있다. 삼각익을 가지는 전투기의 경우 고 받음각으로 비행시 앞전에서 와류가 발생한다. 앞전에서 발생한 와류가 날개의 뒷전까지 흐르지 못하고 박리가 일어나며, 발생한 박리는 와류 붕괴 현상으로 이어진다. 이런 와류 붕괴 현상으로 날개 위에 극심한 비정상 유동이 흐르게 되고, 비정상 유동으로 인해 비행체의 날개가 심하게 흔들리는 현상이 발생한다(Gordnier and Visbal, 2004). 심하게 흔들리는 날개가 다시 날개 위의 비정상 유동을 더욱 심화시켜 고기동이 필요한 전투기의 안정성을 저하시키는 문제를 발생시킨다. 이와 같은 복잡한 물리적 현상을 정확히 예측하기 위해서는 구조와 공력의 정확한 연계 해석이 필요하며, 항공 분야 외에도 가스 터빈, 에어백, 다리 설계, 혈관 해석 등의 다양한 분야에서 공력-구조 연계 해석은 필수적으로 연구가 필요하다.

가장 일반적인 공력-구조 연계해석 방법으로는 유한요소법(Finite Element Method)으로 구조해석을 진행하고, 유한체적법(Finite Volume Method)을 통해 공력해석을 수행하여 두 공간상에서의 해석결과를 경계조건으로 주고는 loosely-coupled 방법이 대표적이다. 공력-구조 연계 해석 문제의 경우 flutter 혹은 대변형 등의 움직임이 많은 문제를 해석하기 때문에 해석 중에 격자의 형상이 크게 변화하는 현상이 발생한다. 이런 경우 정확한 해석을 위해서는 격자를 재생성해야 하며 격자를 재생성에 많은 시간과 노력을 할애해야 한다. 본 연구에서는 이러한 격자의 재생성시 발생하는 어려움을 완화하기 위해서 무격자 기반의 구조해석 기법을 개발하였다.

무격자 해석 기법은 SPH(smooth particle hydrodynamics), RKPM(reproducing kernel particle method), 이동 최소제곱법 (moving least squares approximation), 갤러킨 기법(Galerkin method) 등 여러 가지 무격자 기법들이 있다(Nguyen et al., 2008). 최근에 Katz에 의해 FVM 기반의 공간이산화 대신 최소제곱기반의 무격자해석 기법으로 성공한바가 있다(Katz and Jameson, 2009). 하지만 기하학적보존을 만족하지 않아 물성치 비보전성에 의해 정확한 해석이 불가능하였다. 이에 대한 문제를 Huh가 라그랑지 기법을 도입 및 수학적으로 폐공간을 구성할 수 있는 기하학적 보존 제약 조건을 통해 수치해의 정확성을 상당 부분 개선하였다(Huh et al., 2018). 이는 무격자 기법으로 FVM 기반의 구조해석을 대체할 수 있음을 의미하며, 유동과 구조문제에서 중요하게 여겨지는 격자 생성 및 변형 문제에 상당한 강점이 있다.

본 연구는 코시 모멘텀 방정식(Cauchy’s Momentum equation을 GC-LSM 기법으로 공간차분 하였으며, 시간 전진항은 New Mark Beta(Newmark, 1959) 방법을 활용하였다. 1차원, 2차원, 3차원 문제에 대해서 해석을 진행했으며, 해석해와 비교하여 개발된 기법을 검증하였다. 또한, 질점의 무질서도에 따른 수치해의 정확성을 살펴보기 위해 랜덤한 격자에서도 해석을 진행하여 결과를 비교 검증하였다. 본 연구는 추후 기 보유중인 무격자 유동코드와 결합하여 공력-구조 연계 프로그램 개발을 위한 구조해석 프로그램이며, 구조해석의 자체 해석보다는 공력-구조 연계 해석에서 장점이 있을 것으로 예상된다.

2. 수치기법

본 장에서는 해석에 사용된 지배방정식 및 지배방정식 구성하는 식들을 소개한 후, GC-LSM 기반의 공간차분 방법과 시간전진기법 및 경계 조건을 설명한다.

2.1 지배방정식

본 연구에서 사용한 Cauchy’s Momentum 방정식은 다음과 같다(Mase, 1999).

(1)

ρ a = b + \nabla ∙ σ_{i j}

여기서, b는 체적력, a는 가속도, 𝜌는 밀도 그리고 $σ_{i j}$ 는 응력텐서이다. 식 (1)에서 체적력은 무시하고 외력과 선형 댐핑을 추가하면 식 (2)와 같이 유도된다.

(2)

ρ \frac{\partial v}{\partial t} = \nabla ∙ σ_{i j} + f - c V

(3)

V = [\begin{matrix} u \\ v \\ w \end{matrix}], σ_{i j} = [\begin{matrix} σ_{x x} & τ_{x y} & τ_{x z} \\ τ_{y z} & σ_{y y} & τ_{y z} \\ τ_{z x} & τ_{z y} & σ_{z z} \end{matrix}]

식 (2)에서 f는 외력, V는 속도, c는 댐핑 계수이다.

2.2 구성방정식

응력 변형률 관계식은 등방성 물체 조건에서 Hook’s Law를 사용하였으며 식 (4)와 같다.

(4)

D = \frac{E (1 - v)}{(1 + v) (1 - 2 v)} [\begin{matrix} 1 & \frac{v}{1 - v} & \frac{v}{1 - v} & 0 & 0 & 0 \\ \frac{v}{1 - v} & 1 & \frac{v}{1 - v} & 0 & 0 & 0 \\ \frac{v}{1 - v} & \frac{v}{1 - v} & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1 - 2 v}{2 (1 - v)} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1 - 2 v}{2 (1 - v)} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1 - 2 v}{2 (1 - v)} \end{matrix}]

E는 Young’s Modulus, ν는 Poisson’s 비이다. 변형과 변위 관계식은 Green -Lagrange tensor를 사용하였으며 식 (5)와 같다.

(5)

[\begin{matrix} ε_{x x} \\ ε_{y y} \\ ε_{z z} \\ γ_{x y} \\ γ_{y z} \\ γ_{z x} \end{matrix}] = [\begin{matrix} \frac{\partial d_{x}}{\partial x} \\ \frac{\partial d_{y}}{\partial y} \\ \frac{\partial d_{z}}{\partial z} \\ \frac{1}{2} [\frac{\partial d_{x}}{\partial y} + [\frac{\partial d_{y}}{\partial x}]] \\ \frac{1}{2} [\frac{\partial d_{y}}{\partial z} + [\frac{\partial d_{z}}{\partial y}]] \\ \frac{1}{2} [\frac{\partial d_{z}}{\partial x} + [\frac{\partial d_{x}}{\partial z}]] \end{matrix}]

응력 벡터를 $σ^{T} = [σ_{x x} σ_{y y} σ_{z z} τ_{x y} τ_{y z} τ_{z x}]$ 라고 하고, 변형 벡터를 $ε^{T} = [ε_{x x} ε_{y y} ε_{z z} τ_{x y} τ_{y z} τ_{z x}]$ 라 하면, 식 (3)의 응력 텐서는 다음과 같이 표기할 수 있다. 식 (5)에서 d_x, d_y, d_z는 각 x,y,z 방향의 변위이다.

(6)

σ = D ε

Hook’s Law와 Green Lagrange tensor를 통해 응력을 변위에 대한 2차 미분항으로 표현하였으며, 변위를 미지수로 하는 연립방정식을 구성할 수 있다. 미분 값은 GC-LSM기법을 통해 구하였다.

2.3 무격자 기법

본 연구에서는 Huh가 개발한 GC-LSM(Huh et al., 2018)을 활용하여 식 (6)을 차분한다. 무격자 기법은 질점의 주변 커넥티비티를 이용하여 편미분항을 식 (7)과 같이 계산한다.

(7)

{\frac{\partial ϕ}{\partial x}|}_{i} \approx \sum_{j} a_{i j} (ϕ_{j} - ϕ_{i}) {\frac{\partial ϕ}{\partial y}|}_{i} \approx \sum_{j} b_{i j} (ϕ_{j} - ϕ_{i}) {\frac{\partial ϕ}{\partial z}|}_{i} \approx \sum_{j} c_{i j} (ϕ_{j} - ϕ_{i})

여기서 j는 주변 질점들의 index를 뜻하며, i는 미분값을 근사하고자 하는 점이다. $a_{i j}, b_{i j}, c_{i j}$ 는 무격자 계수로 최소제곱법을 통해 계산한다. 임의의 함수 𝜙에 대해서 Taylor 전개를 통해 질점 $X_{i} = (x_{i}, y_{i}, z_{i})$ 에서의 함수 값을 계산하면 식 (8)과 같다.

(8)

ϕ_{j} = ϕ_{i} + ∆ x_{i j} \frac{\partial ϕ (X_{i})}{\partial x} + ∆ y_{i j} \frac{\partial ϕ (X_{i})}{\partial y} + ∆ z_{i j} \frac{\partial ϕ (X_{i})}{\partial z} + 0 (∆^{2})

식 (8)을 최소화하는 값을 이용해 최소 제곱법으로 무격자 계수를 구한다.

(9)

F = m \in i m i z e \sum_{j = 1}^{n} w_{i j} [∆ ϕ_{i j} - ∆ x_{i j} \frac{\partial ϕ (x_{i})}{\partial x} - ∆ y_{i j} \frac{\partial ϕ (y_{i})}{\partial y} - ∆ z_{i j} \frac{\partial ϕ (z_{i})}{\partial z}]^{2}

식 (9)에서 $w_{i j}$ 는 역거리 가중함수(Mavriplis, 2003)로 식 (10)과 같다.

(10)

w_{i j} = \frac{1}{(∆ x_{i j}^{2} + ∆ y_{i j}^{2} + ∆ z_{i j}^{2})^{1 / 2}}

식 (9)를 통해 구해지는 무격자 계수는 Fig. 1과 같이 볼륨을 구성한다고 했을 때, 닫히는 효과가 나지 않아 물리량 보존이 안되는 현상이 발생한다. 여기서 Geometric Conservation과 1^st order consistency 조건을 추가로 만족하기 위해 Lagrange Multiplier를 사용한다.

(11)

G = \sum_{j}^{n} a_{i j} = 0, H = \sum_{j}^{n} b_{i j} = 0, I = \sum_{j}^{n} c_{i j} = 0

(12)

M_{11} = \sum_{j}^{n} a_{i j} ∆ x_{i j} = 1, M_{12} = \sum_{j}^{n} b_{i j} ∆ x_{i j} = 0, M_{13} = \sum_{j}^{n} c_{i j} ∆ x_{i j} = 0 M_{21} = \sum_{j}^{n} a_{i j} ∆ y_{i j} = 0, M_{22} = \sum_{j}^{n} b_{i j} ∆ y_{i j} = 0, M_{23} = \sum_{j}^{n} c_{i j} ∆ y_{i j} = 0 M_{31} = \sum_{j}^{n} a_{i j} ∆ z_{i j} = 1, M_{32} = \sum_{j}^{n} b_{i j} ∆ z_{i j} = 0, M_{33} = \sum_{j}^{n} c_{i j} ∆ z_{i j} = 0

식 (11)은 Geometric Conservation을 만족하기 위한 식이고 식 (12)은 1^st order consistency를 만족하기 위한 조건이다. 식 (9), (11), (12)을 통해 AX=B 행렬(Huh et al., 2018)을 구성하여 식 (7)에 대한 $a_{i j}, b_{i j}, c_{i j}$ 무격자 계수를 구할 수 있다. 이렇게 구한 무격자 계수는 Fig. 1의 오른쪽 그림처럼 공간이 닫히도록 계수가 구해진다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-01/N0040360108/images/Figure_jcoseik_36_01_08_F1.jpg

Fig. 1.

Schematic of how LSM and GC-LSM construct volume

Fig. 2와 같이 cell interface에서 주변 커넥티비로부터 1차미분항을 구하고, interface에서 구한 1차 미분항을 한번 더 미분하여 $\vec{d} = [d_{x} d_{y} d_{z}]$ 에 대하여 2차분에 대한 식을 구성한다. $\vec{d} = [d_{x} d_{y} d_{z}]$ 는 변위 벡터이다. 1차 미분을 구할 때 커넥티비티의 구성은 i점이 가지고 있는 커넥티비티와 j점이 가지고 있는 커넥티비의 합으로 구성하였으며, interface의 위치는 i점과 j점의 중점으로 하였다. 모든 점에 대해서 위와 같이 식을 구성할 수 있으며, 구성한 식을 행렬 KU=F와 같이 나타낼 수 있다. 강성 행렬 K는 FEM에서 구하는 강성행렬과 같은 역할을 하며 오로지 질점의 분포에 따라 변화한다. 강성행렬 K의 역행렬 K^-1은 PARDISO sparse matrix solver를(MKL, 2020) 사용해 역행렬을 구하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-01/N0040360108/images/Figure_jcoseik_36_01_08_F2.jpg

Fig. 2.

Discretization at cell interface and cell center

2.4 시간전진 기법

Newmark Beta Time Integration을 사용하였으며, 식 (13), (14)에서 β=1/4, γ=1/2를 사용하였다(Newmark, 1959).

(13)

\dot{d_{n + 1}} = {\dot{d}}_{n} + (1 - γ) ∆ t {\ddot{d}}_{n} + γ ∆ t \ddot{d_{n + 1}}

(14)

d_{n + 1} = d_{n} + ∆ t {\dot{d}}_{n} + \frac{∆ t^{2}}{2} ((1 - 2 β) {\ddot{d}}_{n} + d β \ddot{d_{n + 1}})

식 (13), (14)에서 n은 time step을 나타내며, 식 (13), (14)를 식 (2)에 넣고 n step의 값은 F항에 n+1 step의 값은 K항에 넣어 동적 해석을 수행하였다.

2.5 경계조건

동적해석시 고스트 포인트 변위의 예측이 부정확하면 수치 점성에 의해 정확한 수치결과를 얻을 수 없다. 그러므로 고정 벽면과 자유 벽면에서 고스트 포인트의 위치를 부여하는 것은 구조해석에 있어 매우 중요한 일이다. 본 연구에서는 아래와 같은 방법으로 경계조건을 부여하였다.

(16)

d = d

(17)

σ ∙ n = \bar{t}

자유 벽면에서는 Fig. 2의 오른쪽 그림과 같이 벽면 한쪽에 힘이 들어오는 경우 1차 미분 값이 아닌 외력 값을 주었다. Fig. 3 왼쪽 그림은 고정 벽면에서 변위가 0이 나올 수 있도록 고스트 포인트에는 반대 방향의 변위를 주었다. 그리고 자유 벽면에서는 벽면 첫 번째 점에서의 미분 값으로 외삽을 통해 고스트점의 변위를 계산하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-01/N0040360108/images/Figure_jcoseik_36_01_08_F3.jpg

Fig. 3.

Boundary condition at fixed and traction wall

2.4 수치해석 결과 및 고찰

2.4.1 1차원 Bar 해석

Fig. 4는 길이가 L인 bar가 길이에 비례해서 힘을 받는 상태이다. 왼쪽 벽면에서는 변위가 0이며, 오른쪽은 자유단이다. 이 문제의 경우 지배 방정식 및 경계조건은 다음과 같이 주어진다(Dolbow and Belytschko, 1998).

(18)

E u_{x x} + x = 0

(19)

u (0) = 0

(20)

u_{x} (1) = 0

식 (18), (19), (20)에 대한 해석해는 식 (21)과 같이 주어진다.

(21)

u (x) = \frac{1}{E} [\frac{1}{2} x - \frac{x^{3}}{6}]

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-01/N0040360108/images/Figure_jcoseik_36_01_08_F4.jpg

Fig. 4.

A bar of unit length subjected to a linear body force and fixed at wall

Fig. 5에서 볼 수 있듯이 변위 값과 변형률이 해석해와 비교했을 때 정확히 따라가는 것을 볼 수 있다. 또한 Table 1과 같이 격자 수를 늘리면서 에러를 비교하여 공간 정확도를 계산하였다. 에러는 식 (22)와 같이 구하였다.

(22)

L_{2} (e n e r g y) = \sqrt{\frac{1}{2} \sum E (u_{x, ν m} - u_{x, e x a c t})^{2}}

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-01/N0040360108/images/Figure_jcoseik_36_01_08_F5.jpg

Fig. 5.

Comparison of Meshless Method and exact solution at 41node

Table 1.

Grid Number for Grid Convergence Test

Grid Num

121

161

201

251

격자에 따른 error norm을 로그로 표기하여 공간 정확도를 구하였다. Fig. 6에서 격자가 조밀해짐에 따라 무격자 해석 기법이 정확히 2차 정확도를 갖는 것을 볼 수 있으며, 이는 모멘텀 방정식을 해석할 수 있음을 보여준다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-01/N0040360108/images/Figure_jcoseik_36_01_08_F6.jpg

Fig. 6.

Gird Convergence Test on 1-D Bar Problem

2.4.2 2차원 Low Aspect Ratio Beam해석

끝단에서 전단응력을 받고 있는 Timoshenko Beam해석을 수행하였다. 해석개요는 Table 2와 같다.

Table 2.

Parameters of Timoshenko Beam

Parameters
Length	48
Height	12
Young’ Modulus	3e7
Poisson’s Ratio	0.3
Density	1
P	1000

Timoshenko Beam의 해석해는 다음과 같다(Timoshenko and Goodier, 1970).

(23)

u_{x} = - \frac{P y}{6 E 1} [(6 L - 3 x) x + (2 + v) (y^{2} - \frac{D^{2}}{4})]

(24)

u_{y} = \frac{P y}{6 E 1} [3 v y^{2} (L - x) x + (4 + 5 v) \frac{D^{2} x}{4} (3 L - x) x^{2}]

식 (23), (24)의 해석해와 무격자 기법으로 해석한 결과를 Table 3에 비교하였다. Table에서 확인할 수 있듯이 해석해를 잘 따라가며, 격자가 조밀해짐에 따라 높은 정확성을 보여준다. 또한 Fig. 7에서 contour 양상도 물리적으로 잘 나오는 것을 확인하였다.

Table 3.

Comparison of Tip Displacement on Short Beam

Num Points	exact	Meshless	error(%)
36	-0.0089	-7.4905E-3	15.84
64	-0.0089	-8.3497E-03	6.18
100	-0.0089	-8.6816E-03	2.45
144	-0.0089	-8.8198E-03	0.90
196	-0.0089	-8.8812E-03	0.21
256	-0.0089	-8.8898E-03	0.11

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-01/N0040360108/images/Figure_jcoseik_36_01_08_F7.jpg

Fig. 7.

Contour of 2D Beam

2.4.3 2차원 Dynamic Problem

끝단에 힘을 주어 static문제를 해석한 후에 힘을 풀어 자유 진동하는 문제를 해석하였다. Static문제를 검증 후에 시간에 따른 tip의 변위를 이론 값과 비교하였다. 해석 개요는 Table 4와 같으며, 해석해와의 비교를 위해 Poisson’s Ratio는 0으로 설정했다. 본 문제에서stactic과 dynamic의 이론값은 다음과 같다(Xia et al., 2007).

(25)

u_{y} = \frac{4 F l^{3}}{E d h^{3}} = - 0.1

(26)

\frac{y (t)}{y_{s s}} = e^{- ξ ω_{n} t} [\cos ω_{n} t + (\frac{ξ}{\sqrt{1 - ξ^{2}}}) \sin ω_{d} t]

Table 4.

Parameters of Timoshenko Beam

Parameters
Length	20
Height	2
Young’ Modulus	10e6
Poisson’s Ratio	0
Density	2600
P	250
△t	0.2

끝단에서 최대 변위 값은 0.1이며 이론 주파수는 0.05Hz, 이론 주기는 20s가 된다. Table 5에 무격자 기법 해석 결과와 해석해를 비교하였다.

Table 5.

Comparison of Tip Displacement on Long Beam

Num Points	exact	Meshless	error(%)
40	-0.1	1.0875E-03	101.09
90	-0.1	-8.0525E-02	19.48
160	-0.1	-9.2451E-02	7.55
250	-0.1	-9.7198E-02	2.80
360	-0.1	-9.9203E-02	0.80
810	-0.1	-9.9944E-02	0.06

Fig. 8은 끝단에서의 변위를 해석해와 비교한 것이다. 해석해와 비교했을 때 약 200초 동안 lagging 없이 이론해를 따라가는 것을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-01/N0040360108/images/Figure_jcoseik_36_01_08_F8.jpg

Fig. 8.

Comparison of 2D Dynamic Free Vibration Problem

이는 적분없이 2차미분을 구하기 때문에 2^nd order 정확도를 가지고 있어 수치 점성이 낮기 때문에 나올 수 있는 결과라고 생각한다.

2.4.4 3차원 Random Point Static Problem

무격자성을 확인하기 위해 비정렬 질점계에서 해석을 수행하였다. 2.4.2에서 해석한 2차원 문제를 3차원으로 확장하여 점들에 랜덤 계수를 곱하여 점을 흩으려 해석을 수행하였다. 해석 개요는 Table 2와 같으며 다음과 같은 식으로 격자에 랜덤성을 부여하였다.

(27)

x_{r a n d o m} = x_{r e g u l a r} + κ ∙ d_{x} ∙ R N_{x}

(28)

y_{r a n d o m} = y_{r e g u l a r} + κ ∙ d_{y} ∙ R N_{y}

(29)

z_{r a n d o m} = z_{r e g u l a r} + κ ∙ d_{z} ∙ R N_{z}

여기서, 𝜅는 랜덤성을 나타내는 계수이고 0.2를 주었다. d는 격자 간격, RN은 -0.5~0.5 사이의 랜덤 수이다.

Table 6과 Table 7을 비교하여 보면 질점계가 흐트러질 경우 비슷한 정확도를 따라가기 위해서 더욱 많은 점을 필요로 하게 되지만, 질점이 많을수록 이론값에 다가가는 것을 볼 수 있다. 또한 Fig. 9와 Fig. 10에서 변위 contour에서도 확인할 수 있듯이 질점의 배치가 임의로 있어도 해석하는데 문제가 없음을 확인하였다. 본 결과를 통해서 GC-LSM 기반의 무격자 기법의 성공을 확인할 수 있었다.

Table 6.

Comparison of Tip Displacement 3D Structured Beam

Num Points	exact	Meshless	error(%)
32	-0.00899	-5.4838E-03	38.38%
256	-0.0089	-8.1707E-03	8.19%
864	-0.0089	-8.7401E-03	1.80%
2048	-0.0089	-8.8467E-03	0.60%
4000	-0.0089	-8.8742E-03	0.29%
6912	-0.0089	-8.8829E-03	0.19%

Table 7.

Comparison of Tip Displacement 3D Unstructured Beam

Num Points	exact	Meshless	error(%)
864	-0.0089	-6.8067E-03	23.52%
2048	-0.0089	-8.1805E-03	8.08%
4000	-0.0089	-8.2270E-03	7.56%
6912	-0.0089	-8.4348E-03	5.23%
10192	-0.0089	-8.9317E-03	0.36%

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-01/N0040360108/images/Figure_jcoseik_36_01_08_F9.jpg

Fig. 9.

Contour of Structured 3D Beam

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-01/N0040360108/images/Figure_jcoseik_36_01_08_F10.jpg

Fig. 10.

Contour of Unstructured 3D Beam

2.4.5 3차원 Random Point Dynamic Problem

마찬가지로 2.4.3에서 해석한 2차원 문제를 3차원으로 확장하여 정렬 질점계와 비정렬 질점계를 해석하였다. 해석개요는 Table 4와 같으며, 마찬가지로 랜덤 계수𝜅는 0.2를 주었다. Fig. 11은 정령 질점계에서 동적해석 결과이다. 수치점성이 2차원보다는 조금 더 있어 약간의 lagging과 clipping이 발생하는 것을 확인할 수 있지만, 질점 개수가 많아지면 시간 정확도와 공간 정확도가 올라가는 것을 볼 수 있다. Fig. 12는 랜덤 격자에서의 결과이다. 격자의 무질서로 인해 수치 점성이 매우 큰 것을 볼 수 있지만 마찬가지로 질점 개수를 늘렸을 때 해석해를 어렵지 않게 따라가는 것을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-01/N0040360108/images/Figure_jcoseik_36_01_08_F11.jpg

Fig. 11.

Tip Displacement of Structured 3D Beam

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-01/N0040360108/images/Figure_jcoseik_36_01_08_F12.jpg

Fig. 12.

Tip Displacement of Unstructured 3D Beam

3. 결 론

본 연구에서는 GC-LSM 기법을 활용하여 Cauchy’s Momentum 방정식을 차분하였다. 시간전진 기법은 Newmark Beta 기법을 사용했으며, Constitutive Equation으로는 Hook’s Law와 Green- Lagrange tensor를 사용하였다. 위의 식들로 응력을 변위에 대한 2차 편미분함수로 표현하였고, 적분 과정없이 strong form으로 차분하여 해석을 수행했다. 본 연구 방법이 이론적으로 2^nd order 정확도 가지는 것을 확인하였다. 또한 1차원, 2차원 3차원 문제에 대해서 검증을 수행하였고, 정적 동적 문제에 대해서 이론해와 비교하여 정확하게 해석하는 것을 확인하였다.

본 기법은 기존의 FEM에 비해 질점을 많이 필요로 하고 비정렬 질점계에서 수치점성이 비교적 크지만, 본 연구의 목표는 구조해석에서 끝나는 것이 아닌 기 보유중인 무격자 유동 프로그램과 연계하여 무격자 기반의 구조-공력 연계 해석 프로그램을 만드는 것이므로 추후에 FSI문제 해석시에 많은 장점이 있을 것으로 예상된다. 첫째로는 같은 해석기법을 쓰기 때문에 경계면에서 처리해야 하는 일이 적을 것으로 예상된다. 둘째로는 질점의 재생성 및 재구성이 용이하여 구조물의 대변형 문제에 대해 장점이 있을 것으로 예상된다.

본 연구를 통해서 무격자 기반의 공력-구조 연계해석 프로그램 개발의 기틀을 다졌으며, 추후 연구로는 비선형 해석을 추가로 개발 후에 기 보유 중인 유동 해석코드와 연계할 예정이다.

Acknowledgements

본 연구는 국토교통부의 재원으로 국토교통과학기술진흥원 철도기술연구사업(22RTRP-B146018-05)의 지원으로 작성되었습니다.

References

Dolbow, J., Belytschko, T. (1998) An Introduction to Programming the Meshless Element Free Galerkin Method, Arch. Comput. Methods Eng., 5, pp.27~241. 10.1007/BF02897874

Gordnier, R.E., Visbal, M.R. (2004) Computation of the Aeroelastic Response of a Flexible Delta Wing at High Angles of Attack, J. Fluids & Struct., 19, pp.785~800. 10.1016/j.jfluidstructs.2004.04.008

Huh, J.Y, Rhee, J.S., Kim, K.H., Jung, S.Y. (2018) New Least Squares Method with Geometric Conservation Law (GC-LSM) for Compressible Flow Computation in Meshless Method, Comput. & Fluids, 172, pp.122~146. 10.1016/j.compfluid.2018.06.010

Intel, M.K.L. (2020) Intel Math Kernel Library, Developer Reference.

Katz, A., Jameson, A. (2009) A Comparison of Various Meshless Schemes Within a Uniﬁed Algorithm, 47th AIAA Aerospace Sciences Meeting Including The New Horizons Forum and Aerospace Exposition, Orlando, Florida. 10.2514/6.2009-596

Mase, G.T., Smelser, R.E., Mase, G.E. (1999) Continuum Mechanics for Engineers, CRC press, p.398. 10.1201/9780367803230PMC5926170

Mavriplis, D. J. (2003) Revisiting the Least-Squares Procedure for Gradient Reconstruction on Unstructured Meshes, AIAA 16th Computational Fluid Dynamics Conference, Orlando, FL, AIAA, pp.2003~3986. 10.2514/6.2003-3986

Newmark, N.M. (1959) A Method of Computation for Dynamic, J. Eng. Mech. Div., 85(3), pp.67~94. 10.1061/JMCEA3.0000098

Nguyen, V.P., Rabczuk, T., Bordas, S., Duflot, M. (2008) Meshless Methods: A Rreview and Computer Implementation Aspects, Math & Comput. Simul., 79, pp.763~813. 10.1016/j.matcom.2008.01.003

Slotnick, J.P., Khodadodoust, A., Alonso, J., Darmofal, D., Gropp, W., Lurie, E., Mavriplis, D. J. (2014) CFD Vision 2030 Study: A Path to Revolutionary Computational Aerosciences, No. NF1676L-18332, NASA Technical Report, p.20.

Timoshenko, S.P., Goodier, J.N. (1970) Theory of Elasticity (Third ed.), New York, 23, McGraw Hill.

Xia, G.H., Zhao, Y., Yeo, J.H., Lv, X. (2007) A 3D Implicit Unstructured-Grid finite Volume Method for Structural Dynamics, Comput. Mech., 40, pp.299~312. 10.1007/s00466-006-0100-7

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

Preview

Development of Meshless Method Using Least-Squares Method with Geometric Conservation Law for Structural Dynamic Analysis

ABSTRACT

MAIN

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

Fig. 1.

Schematic of how LSM and GC-LSM construct volume

Fig. 2.

Discretization at cell interface and cell center

(13)

(14)

(16)

(17)

Fig. 3.

Boundary condition at fixed and traction wall

(18)

(19)

(20)

(21)

Fig. 4.

A bar of unit length subjected to a linear body force and fixed at wall

(22)

Fig. 5.

Comparison of Meshless Method and exact solution at 41node

Table 1.

Grid Number for Grid Convergence Test

Fig. 6.

Gird Convergence Test on 1-D Bar Problem

Table 2.

Parameters of Timoshenko Beam

(23)

(24)

Table 3.

Comparison of Tip Displacement on Short Beam

Fig. 7.

Contour of 2D Beam

(25)

(26)

Table 4.

Parameters of Timoshenko Beam

Table 5.

Comparison of Tip Displacement on Long Beam

Fig. 8.

Comparison of 2D Dynamic Free Vibration Problem

(27)

(28)

(29)

Table 6.

Comparison of Tip Displacement 3D Structured Beam

Table 7.

Comparison of Tip Displacement 3D Unstructured Beam

Fig. 9.

Contour of Structured 3D Beam

Fig. 10.

Contour of Unstructured 3D Beam

Fig. 11.

Tip Displacement of Structured 3D Beam

Fig. 12.

Tip Displacement of Unstructured 3D Beam

Acknowledgements

References