FESADev/docs/Paper/유한요소해석법을이용한쉘구조물의동적좌굴해석/유한요소해석법을이용한쉘구조물의동적좌굴해석_001.md

![Figure](images/chunk-001-fig-001.jpg)

#### 저작자표시 - 비영리 - 변경금지 2.0 대한민국

#### 이용자는 아래의 조건을 따르는 경우에 한하여 자유롭게

l 이 저작물을 복제, 배포, 전송, 전시, 공연 및 방송할 수 있습니다.

#### 다음과 같은 조건을 따라야 합니다 :

![Figure](images/chunk-001-fig-006.jpg)

저작자표시. 귀하는 원저작자를 표시하여야 합니다.

![Figure](images/chunk-001-fig-008.jpg)

비영리. 귀하는 이 저작물을 영리 목적으로 이용할 수 없습니다.

![Figure](images/chunk-001-fig-010.jpg)

변경금지. 귀하는 이 저작물을 개작, 변형 또는 가공할 수 없습니다.

- l 귀하는, 이 저작물의 재이용이나 배포의 경우, 이 저작물에 적용된 이용허락조건 을 명확하게 나타내어야 합니다.
- l 저작권자로부터 별도의 허가를 받으면 이러한 조건들은 적용되지 않습니다.

#### 저작권법에 따른 이용자의 권리는 위의 내용에 의하여 영향을 받지 않습니다 .

이것은 이용허락규약 (Legal Code)을 이해하기 쉽게 요약한 것입니다.

Disclaimer

![Figure](images/chunk-001-fig-016.jpg)

![Figure](images/chunk-001-fig-017.jpg)

#### 공학석사학위 청구논문

### 유한요소해석법을 이용한 쉘 구조물의 동적 좌굴 해석

## Dynamic Buckling Analysis of Shell Structures using Finite Element Method

![Figure](images/chunk-001-fig-021.jpg)

인하대학교 대학원 항공우주공학과 이 희 준

#### 공학석사학위 청구논문

### 유한요소해석법을 이용한 쉘 구조물의 동적 좌굴 해석

Dynamic Buckling Analysis of Shell Structures using Finite Element Method

2012 년 2 월

지도교수 조 진 연 이 논문을 석사학위 논문으로 제출함

인하대학교 대학원 항공우주공학과 이 희 준

### 이 논문을 이희준의 석사학위논문으로 인정함

### 2012 년 2 월

![Figure](images/chunk-001-fig-031.jpg)

| 주심 | 김 | 기 | 욱 |  |  |
| --- | --- | --- | --- | --- | --- |
|  |  |  |  |  |  |
| 부심 | 조 | 진 | 연 |  |  |
|  |  |  |  |  |  |
| 위원 | 이 | 승 | 수 |  |  |

# 요 약

매개 변수 공진으로도 알려진 동적 좌굴 현상은 구조물이 축 방향의 동적 압축 하중을 받을 때 발생하는 동적 불안정 현상으로서 구조물에 심각한 파손을 유발할 수 있다. 특히 초음속으로 운동하는 항공기나 탄도 미사일, 발사용 로켓과 지구 대기권 재돌입체 그리고 초공동 수중운동체 와 같이 동적 압축 하중을 받는 구조물을 설계할 때 구조물의 동적 좌굴 거동을 예측하여 설계하는 것이 중요하다. 이에 본 논문에서는 축 방향의 동적 압축 하중을 받는 쉘 구조물에 대해 동적 좌굴 해석을 하기 위한 유한요소해석 프로그램을 개발하였으며, 선형/비선형 정적 해석과 진동 및 정적 좌굴 해석을 통해 프로그램에 사용된 쉘 요소의 신뢰성을 확인 하였다. 또한 다양한 모델에 대한 동적 좌굴 해석 결과를 이론적인 해나 실험을 통해 나온 결과와 비교함으로써 본 프로그램의 타당성을 검증하 였다.

### ABSTRACT

Dynamic buckling, also known as parametric resonance, is one of the dynamic  instability phenomena which may lead to serious failure of structure. It occurs when  compressive dynamic loading of axial direction is applied to the structures.  Therefore it is essential to consider the dynamic buckling behaviors of structures,  especially when the structures is designed to be utilized in compressive dynamic  loading of axial direction such as faster supersonic aircrafts, ballistic missiles,  launcher, re-entry vehicles and supercavitating underwater vehicles. In this study,  the finite element program is developed for dynamic buckling analysis. Linear and  nonlinear static analyses, dynamic analysis and static buckling analysis are  performed to demonstrate the accuracy of the developed program. Also the dynamic  buckling analyses are carried out for various models and the computational results  are verified by comparing with analytical and experimental solutions.

# 목 차

| 요약 | i |
| --- | --- |
| ABSTRACT | ii |
| 목차 | iii |
| 그림 목차 | V |
| 표 목차 | vi |
| 1. 서론 | 1 |
| 2. 이 론 | 3 |
| 2.1. MITC4 Shell Element | 3 |
| 2.2. Geometric Nonlinear Formulation | 7 |
| 2.2.1. Finite Rotation Formulation | 21 |
| 2.2.2. Constitutive Matrix | 22 |
| 2.2.3. Mass Matrix | 25 |
| 2.2.4. 6-DOF Shell Element | 26 |
| 2.3. Buckling Theory | 30 |
| 2.3.1. Static Buckling | 31 |
| 2.3.2. Dynamic Buckling | 32 |
| 2.3.3. Dynamic Buckling Theory of Beam | 34 |
| 3. Numerical Example | 39 |
| 3.1. Linear Static Analysis | 39 |
| 3.1.1. Patch Test | 39 |
| 3.1.1.1. Constant Curvature Patch Test | 40 |
| 3.1.1.2. Constant Shear Patch Test | 41 |
| 3.1.1.3. Constant Twist Patch Test | 43 |
| 3.1.2. Pinched Cylinder | 44 |
| 3.1.3. Hemispherical Shell | 46 |
| 3.2. Geometric Nonlinear Analysis | 48 |
| 3.3. Static Buckling Analysis | 50 |

| 3.3.1. Rectangular Plate Shell | 50 |
| --- | --- |
| 3.3.2. Cylindrical Shell | 52 |
| 3.3.3. Stiffened Square Plate Shell | 55 |
| 3.4. Dynamic Buckling Analysis | 57 |
| 3.4.1. Dynamic Buckling Analysis of Beam | 57 |
| 3.4.2. Dynamic Buckling Analysis of Plate | 59 |
| 3.4.3. Dynamic Buckling Analysis of Stiffened Plate | 61 |
| 4. 결 론 | 63 |
| 참고문헌 | 64 |
| 부 록 | 65 |

![Figure](images/chunk-001-fig-043.jpg)

#### 그림 목차

| Fig. 1 Four-node shell element ···························································································· 4 |
| --- |
| Fig. 2 Interpolation function for the transverse shear strains ··············································· 6 |
| Fig. 3 An arbitrary surface with global Cartesian coordinate system, natural coordinate |
| system and local covariant coordinate system spanned by i g ································· 10 |
| Fig. 4 Local Cartesian coordinate system ·········································································· 23 |
| Fig. 5 Global Cartesian coordinate system and local coordinate system ··························· 26 |
| Fig. 6 Dynamic Buckling Model of Beam ········································································· 34 |
| Fig. 7 Patch Test Mesh ······································································································· 39 |
| Fig. 8 Constant Curvature Patch Test Model ····································································· 40 |
| Fig. 9 Constant Shear Patch Test Model ············································································ 41 |
| Fig. 10 Constant Twist Patch Test Model ··········································································· 43 |
| Fig. 11 Pinched Cylinder Model ························································································ 44 |
| Fig. 12 Pinched Cylinder 1/8 Model ·················································································· 44 |
| Fig. 13 Comparison of Convergence for Pinched Cylinder with ABAQUS ······················ 45 |
| Fig. 14 Comparison of Linear Static Analysis for Pinched Cylinder with ABAQUS ········ 46 |
| Fig. 15 Hemispherical Shell Model ··················································································· 46 |
| Fig. 16 Comparison of Convergence for Hemispherical Shell with ABAQUS ················· 47 |
| Fig. 17 Comparison of Linear Static Analysis for Hemisphrical Shell with ABAQUS ····· 48 |
| Fig. 18 Beam Model for Geometric Nonlinear Analysis ··················································· 48 |
| Fig. 19 Comparison of Geometric Nonlinear Analysis for Beam with ABAQUS ············· 49 |
| Fig. 20 Geometry Change of Beam According to Loads Increase ····································· 49 |
| Fig. 21 Rectangular Plate Shell Model ·············································································· 50 |
| Fig. 22 Cylindrical Shell Model ························································································· 52 |
| Fig. 23 Stiffened Square Plate Shell Model ······································································· 55 |
| Fig. 24 Dynamic Buckling Analysis Model for Beam ······················································· 57 |
| Fig. 25 Dynamic Instability Region of Beam ···································································· 58 |
| Fig. 26 Dynamic Buckling Model for Plate ·········································································59 |
| Fig. 27 Dynamic Instability Region of Plate ······································································ 60 |
| Fig. 28 Dynamic Buckling Analysis Model for Stiffened Plate ········································· 61 |
| Fig. 29 Comparison of Dynamic Instability Region for Stiffened Plate with Plate ··········· 62 |

#### 표 목차

| Table. 1 Constant Curvature Patch Test Results ································································ 41 |
| --- |
| Table. 2 Constant Shear Patch Test Results ······································································· 42 |
| Table. 3 Constant Twist Patch Test Results ······································································· 43 |
| ·45 Table. 4 Comparison of Linear Static Analysis for Pinched Cylinder with Exact Solution |
| ·47 Table. 5 Comparison of Linear Static Analysis for Hemispherical Shell with Exact Solution |
| Table. 6 Comparison of Eigenvalue for Rectangular Plate Shell with ABAQUS ·············· 51 |
| Table. 7 Comparison of Mode Shape for Rectangular Plate Shell with ABAQUS·············· 51 |
| Table. 8 Comparison of Eigenvalue for Cylindrical Shell with ABAQUS ························ 53 |
| Table. 9 Comparison of Critical Buckling Pressure for Cylindrical Shell with Analytic Solution · 53 |
| Table. 10 Comparison of Mode Shape for Cylindrical Shell with ABAQUS ···················· 54 |
| Table. 11 Comparison of Eigenvalue for Stiffened Square Plate Shell with ABAQUS ····· 56 |
| Table. 12 Comparison of Mode Shape for Stiffened Square Plate Shell with ABAQUS ·· 56 |

# 1. 서 론

좌굴(buckling, 挫屈)이란 주로 길이가 그 횡단면의 치수에 비해 큰 구 조물의 양단에 압축하중이 가해졌을 경우, 하중이 어느 크기에 이르면 기 둥이 갑자기 휘는 현상을 말한다. 특히 긴 기둥이나 쉘을 많이 사용하는 항공기, 차량, 선박, 건축물 등의 설계에서 좌굴 문제가 중요하며, 원통형 쉘의 경우 가스와 같은 액체를 저장하는 압력용기로 사용될 뿐만 아니라 잠수함이나 항공기와 같이 외압을 받는 구조물에 사용된다. 이와 같이 다 양한 분야에서 사용되는 쉘 구조물은 그 두께가 길이에 비해 얇기 때문 에 진동이나 좌굴에 취약한 특성을 보인다. 이러한 취약점을 해결하기 위 해서는 우선 구조물에 작용하는 하중의 특성을 파악하고 문제가 발생하 는 부분에 대해 구조물을 어떻게 보강할지 결정해야 한다.

일반적으로 좌굴 문제를 다룰 때 정적 하중만을 고려하지만 특수한 경우에 대해서는 동적 하중도 함께 고려해야만 한다. 예를 들어 초음속으 로 운동하는 항공기나 탄도 미사일, 발사용 로켓과 지구 대기권 재돌입체 그리고 초공동 수중운동체의 경우 빠른 속도로 인해 축 방향으로 매우 큰 동적 압축하중이 작용하게 된다. 이러한 동적 압축하중으로 인하여 발 생하는 좌굴을 동적 좌굴(dynamic buckling) 또는 매개변수 공진(parametric  resonance)이라고 한다. 동적 좌굴이 발생하게 되면 구조물의 횡 방향 운 동이 커지게 되고 이로 인해 구조물의 불안정성이 증가하여 치명적인 손 상을 유발할 수 있다. 그러므로 동적 좌굴 현상을 방지하기 위해서는 이 러한 하중이 작용하는 구조물에 대한 해석을 통해 구조물이 불안정해지 는 영역을 파악하고 이를 고려하여 설계하는 것이 중요하다.

하지만 기존 상용유한요소해석 프로그램의 경우 정적 하중에 대한 정 적 좌굴 해석만 가능하며, 동적 하중 또는 정적 하중과 동적 하중이 동시 에 가해지는 경우에 대한 좌굴 해석이 불가능하다. 이에 본 논문에서는 축 방향의 동적 압축 하중을 받는 쉘 구조물에 대해 동적 좌굴 해석을 하기 위한 유한요소해석 프로그램을 개발하였으며, 선형/비선형 정적 해 석과 진동 및 정적 좌굴 해석을 통해 프로그램에 사용된 쉘 요소의 신뢰

성을 확인하였다. 또한 다양한 모델에 대한 동적 좌굴 해석 결과를 이론 적인 해나 실험을 통해 나온 결과와 비교함으로써 본 프로그램의 타당성 을 검증하였다.

![Figure](images/chunk-001-fig-057.jpg)

### 2. 이론

#### 2.1 MITC4 Shell Element

쉘 구조물을 유한요소 모델로 만들기 위해서 다양한 쉘 요소 가운데 Bathe와 Dvorkin에 의해 개발된 MITC4라는 쉘 요소를 선정하였다. MITC4 쉘 요소는 3차원 솔리드 형상으로부터 쉘 형상을 표현하므로 지배방정식의 유한요소 정식화가 다른 쉘 요소에 비해 간단하다. 또한 쉘 이론을 사용하지 않고 3차원 응력, 변형률을 사용하여 표현되며, 임의의 형상에 대한 두꺼운 쉘과 얇은 쉘 모두 적용 가능하다는 장점이 있다. 그리고 대 변형/회전(작은 변형률)과 재료 비선형에 모두 적용 가능하다. [1]

MITC4 쉘 요소 내부의 임의의 점은 고유 좌표계(natural coordinate system)에 대해 정의 할 수 있으며, 위치 벡터는 식(2.1)과 (2.2)같이 나타낼 수 있다.

$$
t
\mathbf{X} = \sum_{I=1}^{4} N_{I}(\xi, \eta) {}^{0}\mathbf{X}_{I} + \frac{\zeta}{2} \sum_{I=1}^{4} t_{I} N_{I}(\xi, \eta) {}^{0}\mathbf{V}_{I}^{n} \tag{2.1}
$$

임의의 시간 t에서

$$
\mathbf{x} = \sum_{I=1}^{4} N_{I}(\xi, \eta)^{t} \mathbf{X}_{I} + \frac{\zeta}{2} \sum_{I=1}^{4} t_{I} N_{I}(\xi, \eta)^{t} \mathbf{V}_{I}^{n} \tag{2.2}
$$

이 때 $N_t(\xi,\eta)$ 는 형상함수이고, ' $\mathbf{X}_t$ 는 시간 t일 때 노드 I의 좌표를 나타내며, 시간 t=0이면 초기 형상에서 노드 I의 좌표를 나타낸다. 그리고 $t_i$ 는 노드 I의 두께이고 ${}^t\mathbf{V}_i^{r}$ 은 시간 t일 때 노드 I의 두께방향의 법선 벡터(normal vector)를 나타내며, 시간 t=0이면 초기 형상에서 노드 I의

두께방향의 법선 벡터를 나타낸다.

임의의 시간 t에서의 MITC4 쉘 요소의 변위는 식(2.3)과 같이 나타낼 수 있다.

$$
{}^{t}\mathbf{u} = \mathbf{x} - \mathbf{X} = \sum_{I=1}^{4} N_{I}(\xi, \eta) \left( {}^{t}\mathbf{X}_{I} - {}^{0}\mathbf{X}_{I} \right) + \frac{\zeta}{2} \sum_{I=1}^{4} t_{I} N_{I}(\xi, \eta) \left( {}^{t}\mathbf{V}_{I}^{n} - {}^{0}\mathbf{V}_{I}^{n} \right)
= \sum_{I=1}^{4} N_{I}(\xi, \eta) {}^{t}\mathbf{u}_{I} + \frac{\zeta}{2} \sum_{I=1}^{4} t_{I} N_{I}(\xi, \eta) \left( {}^{t}\mathbf{V}_{I}^{n} - {}^{0}\mathbf{V}_{I}^{n} \right)
(2.3)
$$

임의의 시간 t에서의 변위는 시간 t일 때의 형상과 초기 형상의 차로부터 구할 수 있으며, 이 때 ${}^{\prime}\mathbf{u}_{I}$ 는 시간 t일 때 노드 I의 변위를 나타낸다. 그리고 이로부터 변위의 증분은 식(2.4)와 같이 나타낼 수 있다.

![Fig. 1 Four-node shell element](images/chunk-001-fig-073.jpg)

$$
\Delta \mathbf{u} = \sum_{I=1}^{4} N_I(\xi, \eta) \Delta \mathbf{u}_I + \frac{\zeta}{2} \sum_{I=1}^{4} t_I N_I(\xi, \eta) \left( -\alpha_I^{\ \ t} \mathbf{V}_I^2 + \beta_I^{\ \ t} \mathbf{V}_I^1 \right) \tag{2.4}
$$

여기서 $\alpha$ 와 $\beta$ 는 각각 $\mathbf{V}^1$ 과 $\mathbf{V}^2$ 방향 벡터의 회전각이고, 이 때 $\mathbf{V}^1$ 과 $\mathbf{V}^2$ 는 $\mathbf{V}^n$ 으로부터 식(2.5)와 같이 구할 수 있다.

$$
\mathbf{V}^{1} = \frac{\mathbf{e}_{2} \times \mathbf{V}_{n}}{\|\mathbf{e}_{2} \times \mathbf{V}_{n}\|}, \quad \mathbf{V}^{2} = \mathbf{V}_{n} \times \mathbf{V}^{1} \tag{2.5}
$$

이 때 $\mathbf{e}_1$ , $\mathbf{e}_2$ , $\mathbf{e}_3$ 는 전역 직교 좌표계(global Cartesian coordinate system)의 기저(basis)이다. 그리고 만약 $\mathbf{e}_2 \times \mathbf{V}_n \approx 0$ 이라면, $\mathbf{V}^1$ 과 $\mathbf{V}^2$ 는 식(2.6)과 같이 나타낼 수 있다.

$$
\mathbf{V}^1 = \mathbf{e}_3, \quad \mathbf{V}^2 = \mathbf{e}_1 \tag{2.6}
$$

하지만 위와 같이 변위를 정의할 경우 일정한 굽힘 모멘트가 가해질때 요소의 모든 점에서 횡 전단 변형률(transverse shear strain)이 영(零)이될 수 없고, 이로 인해 얇은 형상에 대해 요소의 '잠김현상(locking phenomenon)'이 발생하게 된다. 그러므로 연속체 역학의 가정이 Kirchhoff 쉘의 가정을 포함할지라도 유한요소이산화를 통해 이러한 가정을 표현할수가 없다. 이러한 결점을 해결하기 위해 MITC4 쉘 요소에서는 식(2.7)과 같은 횡 전단 변형률에 대한 보간법을 적용하였다.[1]

![Fig. 2 Interpolation function for the transverse shear strains](images/chunk-001-fig-082.jpg)

$$
\widetilde{\varepsilon}_{\xi\zeta} = \frac{1}{2} (1 + \eta) \widetilde{\varepsilon}_{\xi\zeta}^{A} + \frac{1}{2} (1 - \eta) \widetilde{\varepsilon}_{\xi\zeta}^{C} \tag{2.7a}
$$

$$
\widetilde{\varepsilon}_{\eta\zeta} = \frac{1}{2} (1 + \xi) \widetilde{\varepsilon}_{\eta\zeta}^{D} + \frac{1}{2} (1 - \xi) \widetilde{\varepsilon}_{\eta\zeta}^{B} \tag{2.7b}
$$

#### 2.2 Geometric Nonlinear Formulation

좌굴 해석을 하기 위해서는 기하 강성 행렬(geometric stiffness matrix)이 필요하며, 이를 구하기 위해 MITC4 쉘 요소에 대한 비선형 유한요소 정식화 과정을 수행하였다.

현재 형상에 대한 평형방정식은 식(2.8)과 같다.

$$
\nabla_{\mathbf{X}} \cdot \mathbf{\sigma} + \rho \mathbf{f} = \rho \ddot{\mathbf{u}} \tag{2.8}
$$

식(2.8)에 가상 변위에 대한 현재 형상에서의 가상 일 정리를 사용하면, 식 (2.9)와 같이 쓸 수 있다.

$$
\int_{V} \delta \mathbf{u} \cdot (\nabla_{\mathbf{X}} \cdot \mathbf{\sigma} + \rho \mathbf{f}) dV = \int_{V} \delta \mathbf{u} \cdot \rho \ddot{\mathbf{u}} dV \tag{2.9}
$$

식(2.9)를 인덱스를 사용하여 표현하면 식(2.10)과 같다.

$$
\int_{V} \delta u_{i} \cdot \left( \frac{\partial \sigma_{ij}}{\partial x_{j}} + \rho f_{i} \right) dV = \int_{V} \delta u_{i} \cdot \rho \ddot{u}_{i} dV \tag{2.10}
$$

이 때 식(2.11)을 식(2.10)에 대입하면 식(2.12)와 같은 결과를 얻을 수 있다.

$$
\delta u_i \frac{\partial \sigma_{ij}}{\partial x_j} = \frac{\partial \left(\delta u_i \sigma_{ij}\right)}{\partial x_j} - \frac{\partial \delta u_i}{\partial x_j} \sigma_{ij} \tag{2.11}
$$

$$
\int_{V} \frac{\partial \left(\delta u_{i} \sigma_{ij}\right)}{\partial x_{i}} - \frac{\partial \delta u_{i}}{\partial x_{j}} \sigma_{ij} + \delta u_{i} \rho f_{i} dV = \int_{V} \delta u_{i} \cdot \rho \ddot{u}_{i} dV \tag{2.12}
$$

식(2.12)에 가우스의 발산 정리(Gauss' divergence theorem)를 사용하면 식(2.13)과 같이 나타낼 수 있다.

$$
\int_{\partial V} \delta u_i \sigma_{ij} n_j dA + \int_{V} \left( -\frac{\partial \delta u_i}{\partial x_j} \sigma_{ij} + \delta u_i \rho f_i \right) dV = \int_{V} \delta u_i \cdot \rho \ddot{u}_i dV \tag{2.13}
$$

식(2.13)에서 좌측 첫 번째 항의 면 적분 부분은 식(2.14)와 같이 기하학적 경계조건(geometric boundary condition)과 자연적 경계조건(natural boundary condition)으로 나눌 수 있으며, 이 때 식(2.14)의 좌측 첫 번째 항은 기하학적 경계조건에서 가상변위가 영(零)이므로 생략할 수 있다.

$$
\int_{\partial V_g} \delta u_i \sigma_{ij} n_j dA + \int_{\partial V_m} \delta u_i \bar{t}_i dA + \int_V \left( -\frac{\partial \delta u_i}{\partial x_j} \sigma_{ij} + \delta u_i \rho f_i \right) dV = \int_V \delta u_i \cdot \rho i i_i dV \tag{2.14}
$$

$$
\frac{\partial \delta u_i}{\partial x_j} \sigma_{ij} = \frac{1}{2} \left( \frac{\partial \delta u_i}{\partial x_j} + \frac{\partial \delta u_j}{\partial x_i} \right) \sigma_{ij} = \delta \varepsilon_{ij} \sigma_{ij} \tag{2.14}
$$

그리고 식(2.14)에 식(2.15)를 대입하여 정리하면 현재 형상에서의 가상 일에 대한 식을 식(2.16)과 같이 얻을 수 있다.

$$
\int_{\partial V_m} \delta u_i \bar{t}_i dA + \int_V \delta u_i \rho f_i dV = \int_V \delta u_i \cdot \rho \ddot{u}_i dV + \int_V \delta \varepsilon_{ij} \sigma_{ij} dV \tag{2.16}
$$

식(2.16)에서 좌변은 외력에 의한 가상 일이고, 우변은 내력에 의한 가상 일이다.

비선형 해석을 수행하기 위해서는 Total Lagrangian 기법과 Updated Lagrangian 기법이 있으며, 전자의 방법은 문제를 초기 형상에 대해 정의하는 방법이고, 후자의 방법은 현재 형상에 대해 정의하는 방법이다. 본 논문에서는 전자의 방법을 사용하였으며, 식(2.17)을 사용하여 변형된 현재 형상의 응력과 변형률 등을 초기 형상에 대해 정의하였다.

$$
\int_{V} (A)dV = \int_{V_0} (A)\det(\mathbf{F})dV_0 = \int_{V_0} (A)JdV_0 \tag{2.17a}
$$

$$
\int_{\partial V} (A) \mathbf{n} dA = \int_{\partial V_0} (A) J \mathbf{F}^{-T} \widetilde{\mathbf{n}} dA_0 \quad \text{(Nanson's formular)} \tag{2.17b}
$$

$$
\int_{V} \delta \mathbf{\varepsilon} : \mathbf{\sigma} dV = \int_{V_0} \delta \mathbf{F} : \mathbf{P} dV_0 = \int_{V_0} \delta \mathbf{E} : \mathbf{S} dV_0 \tag{2.17c}
$$

이 때 (A)는 임의의 값을 의미하고, F는 변형 구배(deformation gradient), J는 자코비안(Jacobian)으로 체적변화율을 나타낸다. 그리고 $\epsilon$ 과 $\sigma$ 는 현재 형상에서의 미소 변형률(infinitesimal strain)과 Cauchy 응력을 의미하고, E는 Green-Lagrange 변형률로써 초기 형상에서 정의된다. P는 1차 Piola-Kirchhoff 응력을 나타내고, S는 2차 Piola-Kirchhoff 응력을 나타내며, 두 응력 모두 초기 형상에서 정의된다.

식(2.17)을 사용하여 식(2.16)을 초기 형상에 대해 정의해주면 식(2.18)과 같이 나타낼 수 있다.

$$
\int_{V_0} \delta \mathbf{u} \cdot \rho_0 \ddot{\mathbf{u}} dV_0 + \int_{V_0} \delta \mathbf{E} : \mathbf{S} dV_0 = \int_{\partial V_{0...}} \delta \mathbf{u} \cdot \mathbf{F} \mathbf{S} \widetilde{\mathbf{n}} dA_0 + \int_{V_0} \delta \mathbf{u} \rho_0 \mathbf{f} dV_0 \tag{2.18a}
$$

$$
\int_{V_0} \delta \mathbf{u} \cdot \rho_0 \ddot{\mathbf{u}} dV_0 + \int_{V_0} \delta \mathbf{E} : \mathbf{S} dV_0 = \int_{\partial V_{0_m}} \delta \mathbf{u} \cdot \mathbf{F} \mathbf{S} \widetilde{\mathbf{n}} dA_0 + \int_{V_0} \delta \mathbf{u} \rho_0 \mathbf{f} dV_0 \qquad (2.18a)
\Leftrightarrow \int_{V_0} \delta \mathbf{u} \cdot \rho_0 \ddot{\mathbf{u}} dV_0 + \int_{V_0} \delta \mathbf{F} : \mathbf{P} dV_0 = \int_{\partial V_{0_m}} \delta \mathbf{u} \cdot \mathbf{P} \widetilde{\mathbf{n}} dA_0 + \int_{V_0} \delta \mathbf{u} \rho_0 \mathbf{f} dV_0 \qquad (2.18b)
$$

식(2.18b)는 1차 Piola-Kirchhoff 응력으로 표현한 식이고, 식(2.18a)는 2차 Piola-Kirchhoff 응력으로 표현한 식이다. 본 논문에서는 2차 Piola-Kirchhoff 응력을 사용한 식(2.18a)를 사용하였다.

![Fig. 3 An arbitrary surface with global Cartesian coordinate system(x, y, z), natural coordinate system (\xi, \eta, \zeta) and local covariant coordinate system spanned by \mathbf{g}_i](images/chunk-001-fig-117.jpg)

고유 좌표계(natural coordinate system)에서 공변 기저(covariant basis)는 식(2.19)와 같이 나타낼 수 있고, 변위에 대한 기울기(gradient)를 고유 좌표계와 반공변 기저(contravariant basis)를 사용하여 표현하면 식(2.20)과 같이 표현할 수 있다.

$$
\mathbf{g}_i = \frac{\partial \mathbf{x}}{\partial \xi^i} \tag{2.19}
$$

$$
\nabla \otimes \mathbf{u} = \frac{\partial \mathbf{u}}{\partial \xi^i} \otimes \mathbf{g}^i \tag{2.20}
$$

변형률은 식(2.20)을 사용하여 식(2.21)과 같이 나타낼 수 있다.

$$
\mathbf{\varepsilon} = \frac{1}{2} \left[ \nabla \otimes \mathbf{u} + (\nabla \otimes \mathbf{u})^T \right] = \frac{1}{2} \left[ \frac{\partial \mathbf{u}}{\partial \xi^i} \otimes \mathbf{g}^i + \mathbf{g}^i \otimes \frac{\partial \mathbf{u}}{\partial \xi^i} \right] \tag{2.21}
$$

그리고 변형률 텐서(strain tensor)에서 공변 성분(covariant component)을