Design of Three-Dimensional Functionally Graded Composite Structures to Maximize Stiffness and Natural Frequency

Hak-Yong Jang; Jeonghoon Yoo

doi:10.7734/COSEIK.2026.39.1.41

Preview

Research Paper

Journal of the Computational Structural Engineering Institute of Korea. 28 February 2026. 41-48
https://doi.org/10.7734/COSEIK.2026.39.1.41

Design of Three-Dimensional Functionally Graded Composite Structures to Maximize Stiffness and Natural Frequency

강성과 고유진동수 최대화를 위한 3차원 경사기능 복합재 구조 설계

Hak-Yong Jang¹

Jeonghoon Yoo²^*

장 학용¹

유 정훈²^*

¹Graduate Student, Department of Mechanical Engineering, Yonsei University, Seoul, 03722, Korea

²Professor, Department of Mechanical Engineering, Yonsei University, Seoul, 03722, Korea

¹연세대학교 기계공학과 석사과정

²연세대학교 기계공학과 교수

^{*Corresponding Author}

License (open-access, https://creativecommons.org/licenses/by-nc/4.0/):

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

ABSTRACT

This study proposes a design methodology for three-dimensional functionally graded composite structures (FGCSs) by integrating the concept of functionally graded materials (FGMs) with a topology optimization framework. For the design of multiscale composite structures, anisotropic material properties are evaluated using a representative volume element (RVE)–based machine learning module, combined with rotation matrix calculations based on unit quaternions. In addition, a multi-objective optimization is performed with structural stiffness and natural frequency as objective functions to evaluate the optimized structural configurations and their performance under diverse physical environments. The proposed approach enables systematic design of FGCSs with enhanced stiffness and dynamic characteristics.

Keywords

functionally graded composite structure

topology optimization

representative volume element

machine learning

본 논문에서는 경사기능재료(Functionally Graded Material, FGM) 개념과 위상최적설계 방법을 통해 3차원 경사기능복합재 구조(Functionally Graded Composite Structure, FGCS)의 설계 방법을 제안한다. 다중 스케일 복합재 구조 설계 시 이방성 재료 물성치의 계산은 RVE 기반 ML 모듈과 단위 사원수를 활용한 회전 행렬 계산을 통해 수행된다. 또한, 본 연구에서는 다양한 물리적 환경을 설계에 반영하기 위해 구조물의 강성과 함께 고유진동수를 목적함수로 설정한 다중 목적함수 최적설계를 수행하여 구조물의 최적화 형상과 성능을 분석한다.

키워드

경사기능 복합재 구조

위상최적설계

대표체적요소

머신러닝

MAIN

1. 서 론
2. RVE 기반 ML 모듈 제작
2.1 RVE 기반 단위 구조의 재료 물성치 계산
2.2 단위 구조의 재료 물성치 예측 모델 제작
2.3 단위 구조의 회전에 따른 재료 물성치 계산
3. 최적화 문제의 정식화
3.1 위상최적설계 문제의 정식화
3.2 컴플라이언스 최소화 문제
3.3 고유진동수 최대화 문제
3.3 이중 목적함수 위상최적설계 문제
4. 수치 예제
4.1 단일 목적함수 위상최적설계 문제
4.2 이중 목적함수 위상최적설계 문제
5. 결 론

1. 서 론

본 연구에서는 3차원 경사기능 복합재 구조(Functionally graded composite structure, FGCS)의 설계 방법을 제안한다. 이는 위상최적설계를 통해서 수행되며, 위상최적설계란 제약조건을 만족하며 목적함수의 최대화 또는 최소화를 달성하도록 설계 영역 내 재료 밀도 분포를 통해 구조물의 설계를 수행하는 방법이다(Bendsøe and Sigmund, 2003). 최근에는 적층 제조 기술이 발전함에 따라 복합재 구조물의 설계를 위상최적설계를 통해 수행하고자 하는 연구가 활발히 진행되고 있다.

기존 복합재의 설계는 복합재의 섬유 밀도 및 회전 각도로 배치를 최적화하는 방법을 통해 수행되었다. 섬유 배치를 위한 방법은 대표적으로 FMO(Free material optimization), CFAO(Continuous fiber angle optimization), DMO(Discrete material optimization) 방법이 있다(Bruyneel and Fleury, 2002; Ringertz, 1993; Stegmann and Lund, 2005). 최근에는 복합재 다중 스케일 위상최적설계 방법으로 경사기능재료(Functionally graded material, FGM)의 개념을 복합재 구조 설계에 적용하여 거시적 스케일에서의 재료 분포와 미시적 스케일에서의 미시 구조를 동시에 설계하고자 하는 연구가 진행되고 있다. 최근 FGS의 설계 방법으로 균질화 설계법과 SIMP 방법을 사용하여 다중 스케일 위상최적설계를 수행하였지만, 해당 방법은 재료 물성치 계산에 요구되는 계산 비용이 크다는 단점이 있다(Kim et al., 2020). 이에 대한 대안으로 RVE(Representative volume element) 방법과 ML(Machine learning) 방법을 통해 이방성 복합재의 재료 물성치를 효율적으로 근사하는 방법이 제안되었으며, 이를 통해 복합재의 재료 물성치를 적은 계산 비용으로 근사할 수 있다(Kim et al., 2021). 기존 FGS의 설계 방법론은 재료 물성치 계산 시 높은 계산 비용으로 인해 2차원 설계에 대한 연구가 이루어져 왔다. 그러나 2차원 설계의 경우 설계 자유도가 제한적이기에 실제 구조 시스템의 거동을 반영하는 데 한계가 있다.

본 연구에서는 RVE와 ML 방법을 통해 재료 물성치를 근사하는 방법에 대해 소개한다. 단위 사원수를 이용한 3차원 재료 회전 방법에 대한 방법을 다루고, 이를 위상최적설계 프레임워크에 통합하여 3차원 복합재 구조의 최적 형상을 도출한다. 또한, 구조 시스템 설계 시에는 기계적 성능으로 정적 하중에 대한 강성뿐만 아니라 구조물의 동적 특성도 함께 고려되어야 한다. 이에 본 연구에서는 강성과 기본 고유진동수 두 목적함수를 가중치 합산법을 통해 동시에 고려하여 설계를 수행한다(Marler and Arora, 2010). 이에 대한 결과로 가중치 변화에 따른 복합재 구조물의 형상과 성능을 분석하고자 한다.

2. RVE 기반 ML 모듈 제작

2.1 RVE 기반 단위 구조의 재료 물성치 계산

FGCS의 각 요소는 주기적으로 배열되어 있는 단위 구조로 구성되어 있다. 본 연구에서는 단위 구조를 복합재를 정의하는 최소 체적을 의미하는 RVE로 모델링하여, 비균질 재료 특성을 가진 단위 구조를 균질 재료로 정의한다. 본 연구에서 사용한 단위 구조 모델과 고려 가능한 단위 구조 형상을 Fig. 1과 Fig. 2에 나타내었으며, 단위 구조의 형상은 강화재 두께 설계 변수 W₁, W₂, W₃에 의해 결정되어 변수 조합에 따라 9개의 형상이 설계에 반영된다. 단위 구조를 구성하는 각 재료의 물성치는 Table 1과 같다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-01/N0040390105/images/Figure_jcoseik_39_01_05_F1.jpg

Fig. 1.

Unit structures of RVE model

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-01/N0040390105/images/Figure_jcoseik_39_01_05_F2.jpg

Fig. 2.

Possible shapes of the 3D RVE

Table 1.

Material properties for fiber and matrix

Material	$E$ [pa]	𝜈	𝜌 [kg/m³]
Fiber	10	0.3	1
Matrix	10	0.3	0.2

FGCS의 설계에 있어 단위 구조는 무한히 주기적으로 배열되어 있기에 인접한 단위 구조 사이에는 물리적 연속성이 만족되어야 한다. 이에 각 단위 구조가 서로 맞닿는 경계에서 물리적 연속성을 만족시키기 위해 주기 경계 조건을 적용한다. 주기 경계 조건에 의해 RVE의 각 면에서의 변위장은 다음과 같이 표현된다.

(1)

u_{i} = ε ¯_{i k}^{L} y_{k} + u_{i}^{*}

위 식에서 $ε ¯_{i k}^{L}$ 는 RVE에서 발생된 변형률 tensor 성분으로 전산 해석 시 하중 조건의 역할을 수행한다. i, k는 Cartesian 좌표계를 나타내며, $y_{k}$ 는 k번째 축에 대응하는 경계면의 공간 상의 좌표를 나타낸다. $u_{i}^{*}$ 는 주기 함수로 RVE 사이의 주기성을 표현하는 값이지만, 일반적으로 알 수 없는 값이므로 식 (1)을 경계 조건으로 사용할 수 없다. 이를 위해 $u_{i}^{*}$ 를 제거하는 과정이 필요하다. $u_{i}^{*}$ 를 제거하기 위해 식 (1)을 RVE 모델에서 마주보는 두 면을 나누어 다시 정리하면 다음 식 (2)와 같이 표현된다.

(2)

\{\begin{cases} u_{i}^{j +} = {\bar{ε}}_{i k}^{L} y_{k}^{j +} + u_{i}^{*} \\ u_{i}^{j -} = {\bar{ε}}_{i k}^{L} y_{k}^{j -} + u_{i}^{*} \end{cases}

위 식에서 j는 동일한 Cartesian 좌표축을 법선 벡터로 가지는 두 경계면을 나타내며, +와 -는 서로 마주보고 있는 경계면을 의미한다. RVE가 주기적으로 배열되어 있다고 했을 때, 주기 함수 $u_{i}^{*}$ 는 서로 마주보고 있는 두 면에서 같은 값을 가진다. 따라서 위 두 식의 차를 다시 정리하면 다음과 같이 나타낼 수 있다.

(3)

u_{i}^{j +} - u_{i}^{j -} = ε ¯_{i k}^{L} (y_{k}^{j +} - y_{k}^{j -})

여기서, 3차원을 고려하여 $i$ 와 $k$ 를 다시 정리하면 다음 식 (4)와 같이 표현된다.

(4a)

u_{1}^{j +} - u_{1}^{j -} = ε ¯_{11}^{L} (y_{1}^{j +} - y_{1}^{j -}) + ε ¯_{12}^{L} (y_{2}^{j +} - y_{2}^{j -}) + ε ¯_{13}^{L} (y_{3}^{j +} - y_{3}^{j -})

(4b)

u_{2}^{j +} - u_{2}^{j -} = ε ¯_{21}^{L} (y_{1}^{j +} - y_{1}^{j -}) + ε ¯_{22}^{L} (y_{2}^{j +} - y_{2}^{j -}) + ε ¯_{23}^{L} (y_{3}^{j +} - y_{3}^{j -})

(4c)

u_{3}^{j +} - u_{3}^{j -} = ε ¯_{31}^{L} (y_{1}^{j +} - y_{1}^{j -}) + ε ¯_{32}^{L} (y_{2}^{j +} - y_{2}^{j -}) + ε ¯_{33}^{L} (y_{3}^{j +} - y_{3}^{j -})

위 식에서 $y_{k}^{j +} - y_{k}^{j -}$ 의 값은 RVE 모델의 마주보고 있는 변 사이의 거리이며, 이 식을 RVE 주기 경계 조건으로 사용한다. 그러나 식 (3)은 Poisson 효과를 반영하지 못한다. 이에 본 연구에서는 식 (3)에 Poisson 효과를 반영하는 항을 추가하여 다음 식 (5)와 같이 주기 경계 조건을 수정하였다.

(5)

u_{i}^{j +} - u_{i}^{j -} = ε ¯_{i k}^{L} (y_{k}^{j +} - y_{k}^{j -}) - ν ε ¯_{l k}^{L} (y_{k}^{j +} - y_{k}^{j -})

위 식에서 l은 i와 수직인 다른 두 좌표 축을 의미하며, 𝜈는 재료의 Poisson 비이다. 본 연구에서 비균질 재료인 복합재를 균질 재료로 나타내며, 비균질 재료에 저장된 변형 에너지를 U를 균질 재료의 변형 에너지인 U*와 같다. U와 U*는 다음 식 (6), (7)과 같이 정의된다.

(6)

U = \frac{1}{2} \int σ_{i} ε_{j} d V

(7)

U^{*} = \frac{1}{2} {\bar{σ}}_{i} {\bar{ε}}_{j} V_{R V E}

식 (7)에서 ${\bar{σ}}_{i}$ 과 ${\bar{ε}}_{j}$ 는 RVE 내에 발생된 응력과 변형률의 평균값으로 아래의 식 (8), (9)로 표현된다.

(8)

{\bar{σ}}_{i} = \frac{1}{V_{R V E}} \int σ_{i} d V

(9)

{\bar{ε}}_{j} = \frac{1}{V_{R V E}} \int {\bar{ε}}_{j} d V

구성 방정식을 통해 식 (7)를 다시 정리하면 다음과 같이 나타낼 수 있다.

(10)

U^{*} = \frac{1}{2} {\bar{ε}}_{i} {\bar{C}}_{i j} {\bar{ε}}_{j} V_{R V E}

식 (10)에서 변형 에너지 U*는 유한요소해석을 통해 계산되며, 평균 변형률을 하중 조건으로 부여하면 탄성 텐서의 값을 도출할 수 있다. 탄성 텐서 값을 구하기 위한 하중 조건과 그에 따른 수식은 Table 2에 정리하였다.

Table 2.

Effective properties of RVE according to load case

Load case	Effective material properties
${\bar{ε}}_{i} = 1$ , ${\bar{ε}}_{o t h e r s} = 0$	${\bar{C}}_{i i} = \frac{2 U^{*}}{V_{R V E}}, i = 1 ~ 6$
${\bar{ε}}_{i} = 1$ , ${\bar{ε}}_{j} = 1$ , ${\bar{ε}}_{o t h e r s} = 0$	${\bar{C}}_{i j} = \frac{1}{{\bar{ε}}_{i} {\bar{ε}}_{j}} [\frac{2 U^{*}}{V_{R V E}} - ({\bar{C}}_{i i} {\bar{ε}}_{i}^{2} + {\bar{C}}_{j j}^{-} {\bar{ε}}_{j}^{2})]$ $i \neq j, {i, j} = 1 ~ 3$

2.2 단위 구조의 재료 물성치 예측 모델 제작

본 연구에서는 단위 구조의 설계 변수 변화에 따른 재료 물성치 값을 ML 학습을 통해 제작한 ML 모듈로 근사한다. ML 학습을 위한 학습 데이터로 RVE의 설계 변수 조합에 따른 탄성 tensor 데이터를 확보한다. 이는 RVE의 각 설계 변수인 W₁, W₂, W₃를 0부터 1까지 0.05 간격으로 이산화하여 탄성 tensor 데이터를 얻어 총 9861(21 × 21 × 21)개의 탄성 tensor 데이터를 확보하였다.

ML 모듈은 심층 신경망 학습 방법을 통해 제작하였다. 심층 신경망은 두 개 이상의 은닉층으로 이루어진 신경망으로, 이를 이용할 시 입력과 출력 데이터 사이의 관계가 비선형적이거나 복잡하더라도 예측 정확도가 높다는 장점이 있다. 본 연구에서 심층 신경망을 3개의 은닉층과 각 은닉층에서의 노드 개수를 30개로 구성하였다. 또한 신경망의 각 은닉층에서 tangent sigmoid 함수를 활성화 함수로 사용하여 비선형적 관계의 학습이 효과적으로 이루어지도록 하였다.

심층 신경망의 훈련은 각 노드에서의 가중치와 바이어스의 업데이트를 통해 수행되며, 본 연구에서는 신경망 훈련을 위해 MATLAB에서 지원하는 trainbr 내장 함수를 사용하였다. 이를 사용하면 신경망 훈련 시 발생할 수 있는 과적합을 방지할 수 있는 장점이 있다. 또한, 본 연구에서는 과적합 방지를 위해 전체 데이터를 훈련 데이터 70%, 검증 데이터 15%, 테스트 데이터 15%로 분할하여 훈련을 진행하였다. 심층 신경망 훈련 결과로 얻은 ML 모듈의 예측 정확도는 결정 계수(R-square)로 평가하였다. 결정 계수는 다음 식 (11)과 같이 정의된다.

(11)

R^{2} = \frac{\sum_{i = 1}^{n} {({\hat{y}}_{i} - \bar{y})}^{2}}{\sum_{i = 1}^{n} {(y_{i} - \bar{y})}^{2}}

위 식에서 $y_{i}$ 는 관측값, ${\hat{y}}_{i}$ 는 예측값, $\bar{y}$ 는 관측값의 평균을 의미한다. 또한, 결정 계수가 1에 가까울수록 심층 신경망의 예측 정확도가 높음을 의미한다. 본 연구에서는 신경망 학습을 통해 도출한 모든 ML 모듈에서 결정 계수 값이 0.99 이상의 값을 가지도록 하여 높은 예측 정확도를 가지도록 하였다. ML 모듈 제작 결과로서 설계 변수인 W₃의 값을 0.5로 고정한 후 W₁, W₂ 값의 변화에 따른 재료 물성치로서 ${\bar{C}}_{11}$ , ${\bar{C}}_{22}$ , ${\bar{C}}_{33}$ 의 ML 모듈 예측 결과를 Fig. 3에 나타내었다. 해당 결과를 통해 강화재의 두께가 증가함에 따라 재료 물성치도 연속적으로 증가함을 알 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-01/N0040390105/images/Figure_jcoseik_39_01_05_F3.jpg

Fig. 3.

Contour plots of ${\bar{C}}_{11}$ , ${\bar{C}}_{22}$ , ${\bar{C}}_{33}$ using the DNN model according to the variation of design variable W₁, W₂ of the RVE

2.3 단위 구조의 회전에 따른 재료 물성치 계산

복합재 설계 시 더 높은 설계 자유도를 제공하기 위해서는 재료의 회전도 설계에 고려되어야 한다. 본 연구에서는 단위 사원수를 활용하여 회전에 따른 재료 물성치의 계산을 수행한다. 단위 구조의 3차원 회전은 Fig. 4에 나타낸 것과 같이 회전축과 회전각을 통해 수행된다. 회전축 벡터 $\hat{e}$ 는 3개의 설계 변수 ${\hat{e}}_{y 1}, {\hat{e}}_{y 2}, {\hat{e}}_{y 3}$ 를 통해 정의되며, 각 설계 변수의 범위는 0부터 1로 제한하였다. 그리고 단위 사원수 계산 시 회전축 벡터는 단위 벡터로 정의되어야 하기에 다음 식 (12)와 같이 정규화를 수행한다.

(12)

\hat{e} = ({\hat{e}}_{y 1}, {\hat{e}}_{y 2}, {\hat{e}}_{y 3}) {\hat{e}}_{y 1} = \frac{e_{y 1}}{‖ e ‖}, {\hat{e}}_{y 2} = \frac{e_{y 2}}{‖ e ‖}, {\hat{e}}_{y 3} = \frac{e_{y 3}}{‖ e ‖}

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-01/N0040390105/images/Figure_jcoseik_39_01_05_F4.jpg

Fig. 4.

3D Rotation method using unit quaternion

회전축을 단위 벡터와 회전각으로 나타내었을 때, 단위 사원수 q는 사인 함수와 코사인 함수를 통해서 다음 식 (13)과 같이 정의된다.

(13a)

q = (q_{0}, q_{1}, q_{2}, q_{3})

(13b)

q_{0} = \cos (\frac{θ}{2})), q_{1} = {\hat{e}}_{y 1} \sin (\frac{θ}{2})

(13c)

q_{2} = {\hat{e}}_{y 2} \sin (\frac{θ}{2}), q_{3} = {\hat{e}}_{y 2} \sin (\frac{θ}{2})

이를 활용하여 회전 행렬을 다음 식 (14)와 같이 유도할 수 있다.

(14)

R = [\begin{matrix} q_{1}^{2} + q_{0}^{2} - q_{2}^{2} - q_{3}^{2} & 2 (q_{1} q_{2} - q_{0} q_{3}) & 2 (q_{1} q_{3} + q_{0} q_{2}) \\ 2 (q_{1} q_{2} + q_{0} q_{3}) & q_{2}^{2} + q_{0}^{2 - q_{1}^{2}} - q_{3}^{2} & 2 (q_{2} q_{3} - q_{0} q_{1}) \\ 2 (q_{1} q_{3} - q_{0} q_{2}) & 2 (q_{2} q_{3} + q_{0} q_{1}) & q_{3}^{2} + q_{0}^{2 - q_{1}^{2}} - q_{2}^{2} \end{matrix}]

유도된 회전 행렬을 통해 재료 물성치의 회전 변환을 수행하는 변환 행렬 M을 다음 과정의 식 (15)를 통해서 도출할 수 있다.

(15a)

M = [\begin{array}{l} M_{11} & M_{12} \\ M_{21} & M_{22} \end{array}]

(15b)

M_{11} = [\begin{array}{l} R_{11}^{2} & R_{12}^{2} & R_{13}^{2} \\ R_{21}^{2} & R_{22}^{2} & R_{23}^{2} \\ R_{31}^{2} & R_{32}^{2} & R_{33}^{2} \end{array}]

(15c)

M_{12} = [\begin{array}{l} 2 R_{12} R_{13} & 2 R_{11} R_{13} & 2 R_{11} R_{12} \\ 2 R_{22} R_{23} & 2 R_{21} R_{23} & 2 R_{21} R_{22} \\ 2 R_{32} R_{33} & 2 R_{31} R_{33} & 2 R_{31} R_{32} \end{array}]

(15d)

M_{21} = [\begin{array}{l} R_{21} R_{31} & R_{22} R_{32} & R_{23} R_{33} \\ R_{31} R_{11} & R_{32} R_{12} & R_{33} R_{13} \\ R_{11} R_{21} & R_{12} R_{22} & R_{13} R_{23} \end{array}]

(15e)

M_{22} = [\begin{array}{l} R_{22} R_{33} + R_{23} R_{32} & R_{21} R_{33} + R_{23} R_{31} & R_{21} R_{32} + R_{22} R_{31} \\ R_{32} R_{13} + R_{33} R_{12} & R_{31} R_{13} + R_{33} R_{11} & R_{31} R_{12} + R_{32} R_{11} \\ R_{12} R_{23} + R_{13} R_{22} & R_{11} R_{23} + R_{13} R_{21} & R_{11} R_{22} + R_{12} R_{21} \end{array}]

변환 행렬 M를 활용하여 다음의 관계식 (16)을 통해 회전을 고려한 복합재의 재료 물성치의 계산을 수행할 수 있다.

(16)

\hat{C} = [M] [\bar{C}] [M]^{T}

3. 최적화 문제의 정식화

3.1 위상최적설계 문제의 정식화

위상최적설계 과정 수행 시 거시적 스케일에서의 설계 변수 $ϕ_{ρ}$ 와 미시적 스케일에서 강화재의 두께와 배치를 정의하는 $ϕ_{W_{1}}, ϕ_{W_{2}}, ϕ_{W_{3}}, ϕ_{e y 3}, ϕ_{e y 3}, ϕ_{e y 3}, ϕ_{e y 3}, ϕ_{θ}$ 의 업데이트를 매 iteration마다 수행하며, 설계 변수는 다음 과정의 식 (17)을 통해 mapping 된다.

(17)

[\begin{matrix} ϕ_{ρ} \\ ϕ_{W_{i}} \\ ϕ_{{\hat{e}}_{y j}} \\ ϕ_{θ} \end{matrix}] \overset{Helmholtz filter}{\to} [\begin{matrix} {\tilde{ϕ}}_{ρ} \\ {\tilde{ϕ}}_{W_{i}} \\ {\tilde{ϕ}}_{{\hat{e}}_{y j}} \\ {\tilde{ϕ}}_{θ} \end{matrix}] \to [\begin{matrix} ρ = Heav ({\tilde{ϕ}}_{ρ}) \\ W_{i} = 0.5 \times ({\tilde{ϕ}}_{W_{i}} + 1) \\ {\hat{e}}_{y j} = 0.5 \times ({\tilde{ϕ}}_{{\hat{e}}_{y j}} + 1) \\ θ = \frac{90}{2} \times ({\tilde{ϕ}}_{θ} + 1) \end{matrix}]

먼저 Helmholtz 필터 과정이 수행되었으며, 이를 통해 위상최적설계 시 발생 가능한 Checkerboard 패턴과 메시 독립성 문제를 방지할 수 있다(Lazarov and Sigmund, 2011).

(18)

- \nabla \cdot (R_{c}^{2} \nabla \tilde{ϕ}) + \tilde{ϕ} = ϕ

위 식에서 $R_{c}$ 는 Helmholtz 필터의 필터 반경으로 본 연구에서는 메시의 크기와 동일하게 설정하였다. 해당 과정을 통해 모든 설계 변수 𝜙는 -1과 1사이의 값을 가지는 $\tilde{ϕ}$ 로 mapping 된다. 이후 ${\tilde{ϕ}}_{p}$ 는 Heaviside 함수를 통해 0과 1사이의 값으로 매핑 되며, 나머지 설계 변수는 식 (18)의 마지막 과정을 통해 $W_{i}$ , ${\hat{e}}_{y j}$ 의 경우 0과 1 사이의 값, 𝜃의 경우 0과 90 사이의 값으로 매핑된다.

본 연구에서는 밀도법을 기반으로 위상최적설계를 수행한다. 이에 밀도에 따른 재료의 물성치를 Ramp 방법을 통해 근사하였다. Ramp 방법은 고유진동수 문제에서 밀도가 0에 가까운 공극 영역에서 나타나는 국소적인 고유 모드가 발생하는 것을 방지하는 장점이 있다. Ramp 방법을 통해 최종 복합재의 재료 물성치와 밀도는 아래의 식 (19), (20)을 통해 정의된다.

(19)

C = C_{\min} + \frac{ρ}{1 + q (1 - ρ)} (\hat{C} (ρ, W_{i}, {\hat{e}}_{y j}, θ) - C_{\min})

(20)

D = \frac{ρ_{matrix} V_{matrix} + ρ_{fiber} V_{fiber}}{V_{matrix} + V_{fiber}}

두 식에서 $C_{m i n}$ 은 공극 영역의 탄성 텐서 값, $V_{m a t r i x}$ 와 $V_{f i b e r}$ 는 각각 설계 영역 내 기지재와 강화재의 부피이다.

3.2 컴플라이언스 최소화 문제

컴플라이언스 최소화 문제는 다음 식 (21)과 같이 정의된다.

(21)

Min J (ρ, W_{1}, W_{2}, W_{3}, {\hat{e}}_{y 1}, {\hat{e}}_{y 2}, {\hat{e}}_{y 3}, θ, u) = l (u) s.t. a (u, v) = l (v), u \in U, \forall v \in U G_{1} = \frac{\int_{Ω} ρ d x}{\int_{Ω}^{e} 1 d x} - V C_{total} \leq 0 G_{2} = \frac{\int_{Ω}^{e} ρ \times V_{fiber} d x}{\int_{Ω}^{e} 1 d x} - V C_{fiber} \leq 0

목적함수 $J$ 는 컴플라이언스이며, $V_{t o t a l}$ 과 $V C_{f i b e r}$ 는 각각 재료 전체의 부피와 강화재 부피의 제약조건으로 전체 설계 영역에서 전체 부피와 강화재의 부피 분율(volume fraction)은 $V_{t o t a l}$ 과 $V C_{f i b e r}$ 을 초과하지 않는다. 위 식의 풀이를 위한 평형 방정식은 다음 식 (22)와 같다.

(22)

a (u, v) = \int_{Ω} ε (v) : C (ϕ) : ε (u) d x l (v) = \int_{Ω} t \cdot v d x ε (v) = \frac{1}{2} (\nabla v + (\nabla v)^{T})

t는 traction, u와 v는 변위와 가상 변위 벡터이며 다음 식 (23)의 해공간을 만족한다.

(23)

U = \{v \in H^{1} (Ω)^{n} ∣ v = 0 on Γ_{D}\}

3.3 고유진동수 최대화 문제

고유진동수 최대화 문제는 다음 식 (24)와 같이 정의된다.

(24)

Min J (ρ, W_{1}, W_{2}, W_{3}, {\hat{e}}_{y 1}, {\hat{e}}_{y 2}, {\hat{e}}_{y 3}, θ, u_{i}) = \frac{1}{ω_{i}^{2}} s.t. a (u_{i}, v) = b (u_{i}, v), u_{i} \in U, \forall v \in U G_{1} = \frac{\int_{Ω} ρ d x}{\int_{Ω} 1 d x} - V C_{total} \leq 0 G_{2} = \frac{\int_{Ω} ρ \times V_{fiber} d x}{\int_{Ω} 1 d x} - V C_{fiber} \leq 0

위 식에서 $u_{i}$ 는 $i$ 번째 모드의 고유 벡터, $ω_{i}$ 는 고유진동수이다. $a (u_{i}, v)$ 와 $b (u_{i}, v)$ 는 다음 식 (25)와 같이 정의된다.

(25)

a (u_{i}, v) = \int_{Ω} ε (v) : C (ϕ) : ε (u_{i}) d x b (u_{i}, v) = ω_{i}^{2} \int_{Ω} D (ϕ) u_{i} \cdot ε (v) d x

3.3 이중 목적함수 위상최적설계 문제

컴플라이언스와 기본 고유진동수를 동시에 고려한 이중 목적 위상최적설계 문제는 다음 식 (26)과 같이 정식화 된다.

(26)

Min J = w_{1} {\bar{J}}_{1} + w_{2} {\bar{J}}_{2} J_{1} (ρ, W_{1}, W_{2}, W_{3}, {\hat{e}}_{y 1}, {\hat{e}}_{y 2}, {\hat{e}}_{y 3}, θ, u_{1}) = \frac{1}{ω_{1}^{2}} J_{2} (ρ, W_{1}, W_{2}, W_{3}, {\hat{e}}_{y 1}, {\hat{e}}_{y 2}, {\hat{e}}_{y 3}, θ, u) = l (u) {\bar{J}}_{i} = \frac{J_{i}^{k} - J_{i}^{u t o}}{J_{i}^{init} - J_{i}^{u t o}} w_{1} + w_{2} = 1, 1 \geq w_{1}, w_{2} \geq 0

이중 목적 위상최적설계는 각 목적함수에 가중치를 곱한 후 더하여 하나의 목적함수를 최소화하는 방법인 가중치 합산법을 통해 수행하였다. 이때 식 (26d)를 통해 두 목적함수를 정규화를 통해 두 목적함수가 유사한 스케일을 갖도록 하였다. 식 (26d)에서 $k$ 는 최적설계 과정 중 현재 iteration을 의미하고 $J_{i}^{init}$ 과 은 $J_{i}^{u t o}$ 는 각각 첫 번째 iteration의 목적함수 값과 이상 최적해의 목적함수 값이다.

4. 수치 예제

본 연구에서는 3D Cantilever beam의 설계를 수치 예제로 다룬다. 수치 예제 수행을 위한 하중 조건 및 경계 조건은 Fig. 5와 Fig. 6에 나타내었다. Neumann 경계에는 하중 조건을 부여하였으며, 동일한 영역에 집중 질량을 부여하여 고유진동수 최대화 문제에서 재료가 경계 조건에 집중되는 현상을 방지하였다. 수치 예제에서 부피 제약조건은 $V C_{total}$ = 0.4, $V C_{fiber}$ = 0.12로 설정하여 설계 영역 내 전체 재료의 부피가 40%, 복합재 내 강화재 부피가 30%를 만족하도록 하였다. 최적설계를 위한 알고리즘은 Method of Moving Asymptotes(MMA)를 사용하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-01/N0040390105/images/Figure_jcoseik_39_01_05_F5.jpg

Fig. 5.

Boundary condition of cantilever beam

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-01/N0040390105/images/Figure_jcoseik_39_01_05_F6.jpg

Fig. 6.

Design results of the FGCS for fundamental natural frequency maximization

4.1 단일 목적함수 위상최적설계 문제

단일 목적함수를 고려한 FGCS의 최적화 형상 및 성능을 Fig. 6과 Fig. 7에 나타내었다. 복합재의 기지재를 Fig. 6(a)와 Fig. 7(a)에 나타내었으며, 미시적 스케일에서 강화재의 배치에 대한 결과는 Fig. 6(b)와 Fig. 7(b)에 나타내었다. 해당 그림을 통해 설계 영역 내 위치가 변화함에 따라 강화재의 두께가 연속적으로 변화하는 것을 확인할 수 있다. 추가적으로 Fig. 6과 Fig. 7에 나타낸 구조물의 성능 비교를 통해 하나의 목적함수만 고려할 때 다른 기계적 성능은 감소하는 trade-off 되는 관계를 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-01/N0040390105/images/Figure_jcoseik_39_01_05_F7.jpg

Fig. 7.

Design results of the FGCS for compliance minimization

4.2 이중 목적함수 위상최적설계 문제

본 연구에서는 가중치 합산법의 가중치 $w_{1}$ 을 0부터 1까지 0.1간격으로 변화하여 총 11개의 수치 예제를 수행하였다. 이때 Case 1( $w_{1}$ =1)의 경우에는 기본 고유진동수만 고려한 단일 목적 최적화 결과이며, Case 11( $w_{1}$ =0)의 경우에는 구조물의 강성만 고려한 단일 목적 최적화 결과이다. 이중 목적함수 위상최적설계의 결과로 가중치의 변화에 따른 목적함수의 변화 추세를 나타낸 Pareto 최적해를 Fig. 8에 나타내었다. 컴플라이언스 최소화 문제는 하중 전달 경로를 따라 재료 분포를 집중시키는 반면, 고유진동수 최대화는 진동 모드 증가에 효과적으로 기여하도록 강성과 질량 분포가 유도된다. 이에 두 목적함수 사이의 trade-off 관계가 나타나는 것을 Fig. 8과 Table 3을 통해 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-01/N0040390105/images/Figure_jcoseik_39_01_05_F8.jpg

Fig. 8.

Pareto front of bi-objective topology optimization

Table 3.

Performance of FGCS topology optimization

	Weight $w_{1}$	$ω_{1}$ [Hz]	Compliance [J]
Case 1	1.0	0.2290	72.0650
Case 2	0.9	0.2276	71.3865
Case 3	0.8	0.2223	70.8621
Case 4	0.7	0.2202	70.8921
Case 5	0.6	0.2175	70.2908
Case 6	0.5	0.2158	68.9074
Case 7	0.4	0.2015	63.5666
Case 8	0.3	0.1800	59.8785
Case 9	0.2	0.1700	58.2853
Case 10	0.1	0.1640	57.7910
Case 11	0.0	0.1597	57.9410

또한 Pareto 최적해의 결과로 이상 최적해와 가장 가까운 설계안을 탐색할 수 있다. Fig. 8에서 확인할 수 있듯이 Case 7의 결과가 이상 최적해와 가장 가깝게 도출되었으며, 최적화된 구조물의 형상과 성능을 비교하였을 때, Case 7은 두 단일 목적함수 최적화 결과 사이의 성능을 가지는 것을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-01/N0040390105/images/Figure_jcoseik_39_01_05_F9.jpg

Fig. 9.

FGCS results of bi-objective topology optimization

5. 결 론

본 연구에서는 3차원 경사기능복합재 구조(FGCS)의 설계를 위해 복합재 구성에 따른 재료 물성치 예측 모델을 제작하고 이를 위상최적설계 방법에 통합하여 다중 스케일 위상최적설계 방법을 제시하였다. 제안된 방법은 거시적 설계 변수인 밀도를 기반으로 설계 영역 내 재료의 분포를 결정하고 미시적 설계 변수인 강화재의 두께 변수와 회전 변수를 통해 복합재의 재료 물성치를 근사한다. 강화재 두께에 변화에 따른 재료 물성치의 변화는 RVE기반 ML 모듈을 통해 예측하였다. 추가적으로 복합재 구조 설계 자유도를 높이기 위해 단위 사원수를 활용한 회전 기법을 통해 복합재의 회전을 고려하였다.

제안한 설계 기법을 통해 3D Cantilever beam의 설계를 수행하였다. 수치 예제로 컴플라이언스와 최소화와 기본 고유진동수 최대화 두 목적함수를 가중치 합산법을 통해 동시에 고려하였으며, 가중치 변화에 따라 두 목적함수의 변화 추세를 나타내는 파레토 최적해를 탐색하였다. 또한 이를 통해 이상 최적해와 가장 근접한 설계안을 탐색하였다.

본 연구에서 제안한 방법은 적층 제조 방식의 발전에 따라 복합재 구조 설계에 현실적인 기여를 할 수 있을 것으로 사료되며, 해당 방법을 사용함으로써 다양한 물리적 환경을 고려한 3차원 복합재 설계에 대한 연구로 확장될 수 있다.

Acknowledgements

본 연구는 2022년도 한국연구재단(NRF-2022R1F1A1065611) 및 ㈜ KTE의 지원을 받아 수행되었습니다.

References

Bendsøe, M.P., Sigmund, O. (2003) Topology Optimization: Theory, Methods and Applications, Springer-Verlag, Berlin, pp.1~2.

Bruyneel, M., Fleury, C. (2002) Composite Structures Optimization using Sequential Convex Programming, Adv. Eng. Softw., 33, pp. 697~711.

10.1016/S0965-9978(02)00053-4

Kim, C., Lee, J., Yoo, J. (2021) Machine Learning-Combined Topology Optimiation for Functionary Graded Composite Structure Design, Comput. Methods Appl. Mech. Eng., 387, 114158.

10.1016/j.cma.2021.114158

Kim, D., Lee, J., Nomura, T., Ded,e E.M., Yoo, J., Min, S. (2020) Topology Optimization of Functionally Graded Anisotropic Composite Structures using Homogenization Design Method, Comput. Methods Appl. Mech. Eng., 369, 113220.

10.1016/j.cma.2020.113220

Lazarov, B.S., Sigmund, O. (2011) Filters in Topology Optimization based on Helmholtz‐Type Differential, Int. J. Numer. Methods Eng., 86, 765~781.

10.1002/nme.3072

Marler, R.T., Arora, J.S. (2010) The Weighted Sum Method for Multi-Objective Optimization: New Insights, Struct. Multidiscip. Optim., 41, pp.853~862.

10.1007/s00158-009-0460-7

Ringertz, U. (1993) On Finding the Optimal Distribution of Material Properties, Struct. Optim., 5, pp.265~267.

10.1007/BF01743590

Stegmann, J., Lund, E. (2005) Discrete Material Optimization of General Composite Shell Structures, Int. J. Numer. Methods Eng., 62(14), pp.2009~2027.

10.1002/nme.1259

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

Preview

Design of Three-Dimensional Functionally Graded Composite Structures to Maximize Stiffness and Natural Frequency

ABSTRACT

MAIN

Fig. 1.

Unit structures of RVE model

Fig. 2.

Possible shapes of the 3D RVE

Table 1.

Material properties for fiber and matrix

(1)

(2)

(3)

(4a)

(4b)

(4c)

(5)

(6)

(7)

(8)

(9)

(10)

Table 2.

Effective properties of RVE according to load case

(11)

Fig. 3.

Contour plots of C¯11, C¯22, C¯33 using the DNN model according to the variation of design variable W1, W2 of the RVE

(12)

Fig. 4.

3D Rotation method using unit quaternion

(13a)

(13b)

(13c)

(14)

(15a)

(15b)

(15c)

(15d)

(15e)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

(25)

(26)

Fig. 5.

Boundary condition of cantilever beam

Fig. 6.

Design results of the FGCS for fundamental natural frequency maximization

Fig. 7.

Design results of the FGCS for compliance minimization

Fig. 8.

Pareto front of bi-objective topology optimization

Table 3.

Performance of FGCS topology optimization

Fig. 9.

FGCS results of bi-objective topology optimization

Acknowledgements

References

Contour plots of ${\bar{C}}_{11}$ , ${\bar{C}}_{22}$ , ${\bar{C}}_{33}$ using the DNN model according to the variation of design variable W₁, W₂ of the RVE