1. 서 론
최근 다양한 공학분야에서 구조물의 경량화와 고강도화가 요구되고 있다. 이에 따라 개발된 재료로서 질량대비 강도가 높은 복합재료(Composite materials)가 광범위하게 적용되고 있다. 그러나 이러한 전통적 복합재료는 균질성 혼합체로서 복합재료를 구성하는 재료들의 고유의 특성이 저하될 수 있다(Shaker et al., 2008). 또한 일반적으로 적층형태를 가지므로 재료 사이의 급격한 응력의 변화 등에 따라 층간분리가 발생하기도 하며, 층사이의 과대 변형에 의한 균열의 발생 등을 야기하기도 한다.
전통적 복합재료에서 나타나는 이러한 단점을 보완할 수 있는 대안 중 하나는 기능경사재료(functionally graded materials: FGM)의 적용이다. 기능경사재료는 두 가지의 재료가 점진적으로 전이되어 가며 고체화되는 과정으로 생성되며, 이질재료 연결부의 열전도영역의 연속화와 이질재료 접촉면에서의 응력집중현상 경감 등에 유리하고, 터빈부품, 로켓의 노즐, 화학반응조의 관, 전지 셀 등의 다양한 분야에 적용되고 있다.
Fig. 1은 중립면에 대칭으로 기능경사를 가지는 재료를 도식적으로 보여주고 있다. 중립면과 외면에서의 재료는 각기 다른 재료로서 
, 
의 탄성계수를 가지며, 전이영역에서는 목적하는 함수에 따라 재료상수가 변화된다. 이러한 기능경사는 두 가지의 재료의 체적비가 구조의 두께방향을 따라 연속적으로 변화하면서 구성된다. 복합재료의 한 형태인 기능경사재료는 1984년 일본의 우주계획 프로젝트 중 제안되었다. 기능경사재료는 두 개의 이질 재료가 각 재료의 고유의 특성을 그대로 유지하며 조합된 재료로서, 예를 들어 강과 세라믹을 혼용할 경우 강의 인성과 세라믹 재료의 내화성을 양단 경계에서 그대로 유지할 수 있다.
그러나 기능경사재료의 재료특성치는 두께방향으로 이론적 경사를 만족하도록 제작하는 것이 매우 어려우며, 두 이질 재료의 체적비의 불확실성과 공극 등의 존재로 불확실성을 가지는 것을 피할 수 없다.
이런 부분은 기존에 FGM 평판에 대한 연구에서도 일부 고려되었다. Shaker 등(2008)은 FGM 평판의 재료 불확실성이 고유진동수에 미치는 민감도 해석을 FROM과 SORM을 적용하여 추계론적 유한요소해석 영역에서 제시하였고, Ferrante 등(2005)은 열하중을 받는 FGM 평판에 대해 이질재료의 부피비와 공극에서 나타나는 불확실성을 고려한 해석을 수행하였다. 그러나 이들 연구들은 모두 통계학 영역에 속하는 추계론적 해석법인 몬테카를로 해석(Monte Carlo simulation: MCS)에 기본하고 있다. 잘 알려진 바와 같이 MCS는 불확실인수에 대한 난수 발생 및 이를 이용한 다중의 확정론적 해석으로 구성되며, 해석비용이 매우 높을 뿐 아니라, 난수 발생의 정확성 및 효율성에 크게 영향을 받게 된다.
본 연구에서는 축방향을 따라 보의 폭이 변화되는 FGM 보에서 재료상수의 축방향을 따르는 공간적 불확실성을 고려할 수 있는 비통계학적 추계론적 유한요소해석 정식화를 제시한다. 또한 대상 보의 수치적 모델링을 위하여 Hermite 형상함수에 기반한 2절점 보 유한요소를 정식화하여 사용하였다. 제안한 해석법에 의한 결과는 스펙트럼모사법(spectral representation method; Shinozuka et al., 1991)에 의한 MCS해석을 통하여 검증하였다.
2. 본 론
이 장에서는 기능경사재료 및 이를 적용한 보에 대하여 기술하고, 기능경사재료 보의 해석을 위한 2절점 보요소의 구성에 대해 논한다. 본 논문의 주제인 기능경사재료의 추계론적 해석을 위하여 자기상관함수를 이용한 정식화를 제시하고, 몬테카를로 해석과 비교를 통하여 제안 해석법의 타당성을 제시한다.
2.1 기능경사재료 보의 강성행렬
2.1.1 기능경사재료 보의 모델
기능경사재료로 식 (1)과 같이 중립면(y=0)을 기준으로 대칭이고, 상하 방향으로 지수함수적으로 변화하는 탄성계수를 가지는 것으로 가정한다. 중립면에서의 탄성계수는 
의 값을 가지며 높이 
2에서는 재료지수 
에 따라 제2의 재료에 대한 특정의 값을 가질 수 있다.
(1)
대상 기능경사보는 기하학적 변화도 함께 고려하는데 보의 평면 
에서 
-축을 따라 식 (2)와 같이 보의 폭이 변화하는 변단면 형상을 고려하고자 한다.
(2)
형상인수 β와 축방향 좌표축의 원점의 위치에 따라 지수함수적으로 보의 폭이 변화하는 형태이다.
2.1.2 유한요소 정식화
기술한 재료적 특성과 기하학적 특성을 가지는 기능경사재료 보의 유한요소의 형태는 Fig. 2와 같다.
보가 축방향으로 지수함수로 변화하는 폭을 가지고 있으나 (Fig. 2의 점선), 유한요소 모델 시 요소내의 길이가 상대적으로 짧고, 기하학적 변화는 전체 요소망에 의해서 모사되므로 개별 요소에서는 보의 폭이 선형변화(Fig. 1의 실선)하는 것으로 가정한다. 따라서 폭 
는 다음의 식으로 표현된다.
(3)
요소의 변위성분은 절점당 2개의 자유도를 가정하며, 이에 의한 변위장 
는 Hermite 형상함수에 의한 보간으로 나타낸다.
(4)
여기서, 
는 Hermite 형상함수벡터, 
는 절점변위벡터를 나타내며, 각각 다음과 같다.
(5)
보에서의 곡률 
와 모멘트의 관계 그리고 휨강성인수 
를 고려하면 보의 모멘트는 식 (6)과 같이 얻을 수 있다.
(6)
여기서, 휨강성인수 은 다음과 같다.
(7)
보의 변형에너지가 식 (8)과 같으므로
(8)
최소포텐셜에너지 원리(Cook et al., 1989)를 적용하면 기능경사재료 보의 강성행렬은 식 (9)와 같이 얻을 수 있다.
(9)
그러므로 기능경사재료 보의 요소강성행렬은 
방향으로의 탄성계수 변화 효과를 나타내는 상수 가 곱해진 형태를 가지게 됨을 알 수 있으며, 결국에는 보의 중립면에서의 탄성계수인 
의 함수로 결정된다. 보의 폭 변화의 영향은 
로 고려되며, 적분 후의 최종 요소강성행렬은 요소 양단의 폭 
로 표현된다.
2.2 추계론적 해석을 위한 정식화
본 논문에서는 기능경사재료의 중립면에서의 탄성계수 의 축방향을 따르는 공간적 불확실성을 고려한 추계론적 유한요소해석에 대하여 고찰하고자 한다. 즉, 탄성계수 은 보의 축방향좌표 
의 함수가 되며, 높이 방향의 탄성계수 
는 식 (1)에 따라 에 종속된다.
2.2.1 불확실성을 가지는 강성행렬의 구성
식 (9)로 표현된 기능경사재료 보의 강성행렬에서 중립면 탄성계수를 보의 축방향에 대한 임의성을 포함하는 식으로 표현하면 다음과 같다(Noh, 2008).
(10)
는 중립면 탄성계수의 평균치 또는 확정론적 값을 의미하며, 
는 추계장함수(stochastic function)이다. 이를 식 (7)에 대입하면 휨강성인수는 평균항과 변동항으로 구분된다.
(11)
그러므로 식 (9)의 요소강성행렬 
는 평균 탄성계수만을 포함하는 평균강성행렬 
와 추계장이 곱해진 변동강성행렬 
의 합으로 나뉘어진다.
(12)
식 (12)의 변동강성행렬은 추계장함수 의 함수로서 추계적분(stochastic integration)에 의한 확률변수를 포함하고 있음을 알 수 있다. 추계적분은 미지의 함수인 추계장함수 와 기지의 함수 
의 곱에 대한 적분으로 정의되며, 보통 식 (13)과 같이 확률변수 
로 정의된다(Noh, 2011).
(13)
식 (9)로 얻어지는 강성행렬은 기능경사보의 길이 
과 양단에서의 폭 
그리고 휨강성인수 평균치인 
의 함수로 나타나며, 모두 6개의 독립적인 다항식으로 구성된다: 
, 
1, 2. 그러므로 식 (12)의 변동 강성행렬은 개별 유한요소별로 6개의 독립적 확률변수를 가진다. 즉, 
. 그러므로 요소강성행렬 또한 6개의 확률변수를 가지게 된다.
2.2.2 변위의 확률모멘트 산정
구조 거동의 분산도를 산정하기 위해서는 변위에 대한 1차 모멘트(평균)와 2차 모멘트(공분산)의 결정이 요구된다. 이를 위하여 먼저 변위벡터를 평균확률변수에 대해 선형 Taylor 전개하고, 이를 용하여 응답의 1, 2차 모멘트를 구한다.
1) 변위벡터의 선형 Taylor 전개
앞에서 기술한 바와 같이 요소강성행렬이 요소별로 개별적인 6개의 확률변수의 함수로 표현되므로 구조계 강성행렬과 하중벡터와의 연산을 통한 변위벡터 또한 확률변수의 함수로 가정할 수 있다. 따라서 변위벡터에 대한 확률변수의 평균치를 중심으로 한 전개식을 다음과 같이 기술한다.
(14)
여기서, 
는 구조모델 내 유한요소의 수, 
은 개별 요소내의 확률변수의 개수로 6의 값을 가진다. 
의 평형방정식에 대한 편미분을 통하여 확률변수에 대한 변위의 편미분은 구조계 강성행렬에 대한 편미분으로 대치될 수 있으므로, 식 (14)는 식 (15)로 변형된다.
(15)
식 (14), (15)의 아래첨자로 사용된 
는 평균확률변수벡터에서 계산된 값을 나타낸다.
2) 변위의 1차 및 2차 모멘트
변위의 1차 모멘트(평균)는 식 (15)에 대한 평균연산으로 얻을 수 있다. 이 경우 확률변수의 평균이 
이므로 식 (15)의 두 번째 항이 사라지며, 따라서 평균 변위벡터는 다음과 같이 구해진다. 여기서 
은 평균연자를 의미한다.
(16)
즉, 평균변위벡터는 평균강성행렬과 하중벡터 사이의 연산에 의해 얻어지는 결정론적 해석의 결과와 동일한 값을 나타낸다.
변위에 대한 2차모멘트는 공분산연산을 통하여 얻을 수 있다. 즉,
(17)
식 (15)와 (16), 그리고 확률변수의 평균이 0이라는 조건을 이용하면, 식 (17)은 다음과 같이 정리된다.
(18)
식 (18)에서 구조계 강성행렬의 확률변수에 대한 편미분과 확률변수의 곱에 대한 합은 요소 변동강성행렬로 귀결되며, 이를 통하여 식 (18)은 다음과 같이 쓸 수 있다.
(19)
여기서, 
이며, 식 (19) 내의 평균연산은 다음과 같다.
(20)
여기서, 
. 그러므로 식 (16), (19), (20)을 이용하여 변위에 대한 1차 및 2차 모멘트를 구하고, 구조응답의 분산계수(coefficient of variation: COV)를 아래와 같이 산정한다.
(21)
분산 
는 식 (19)에 주어진 공분산행렬에서 해당 자유도의 대각요소에 해당한다.
2.3 몬테카를로 해석
기능경사재료 보에 대해 제시한 해석법은 해석 결과의 일관성을 통해 그 자체로 타당성을 보일 수도 있겠으나, 추계론적(또는 확률론적)해석의 근간인 몬테카를로 해석을 통한 검증이 보다 일반적이다. 본 연구에서는 이를 위하여 스펙트럼모사법에 기초하여 랜덤 프로세스를 생성하였다. 스펙트럼모사법에 의한 랜덤 프로세스 생성은 다음의 식에 따른다.
(22)
여기서, 
,
,
이며, 
은 [0, 2
]내에 등분포하는 임의위상각(random phase angle)이다. 생성되는 랜덤 프로세스는 N이 증가할수록 정규분포로 수렴하며, 본 연구에서는 N=512를 사용하였다.
2.4 수치예제
본 논문에서는 길이 
20.0m, 높이 
0.3m, 초기 보폭 
1.0m의 기능경사재료 단순보에 대한 해석 결과를 제시하고자 한다. 작용하중은 두 경우를 고려한다. 첫 번째는 중앙점 집중하중 P=10.0kN(하항)(기호, P)이고, 두 번째는 중앙점 하중에 더하여 M=250.0kNm(반시계방향)의 모멘트하중이 추가로 작용하는 경우(기호, P&M)이다.
몬테카를로 해석에서는 모든 경우에 10,000개의 샘플을 생성/사용하며, 해석시간 감소 및 정확도의 개선을 위해 국부평균법(local average scheme)을 적용한다. 제안한 정식화에 의한 수치해석 및 몬테카를로 해석에서는 식 (23)의 자기상관함수(auto-correlation function)와 스펙트럼밀도함수(spectral density function)를 사용하였다.
(23)
여기서, 
는 추계장의 분산, 
는 상관관계거리, 
는 영역 내 두 점사이의 상대거리를 나타낸다.
Figure 3
Examples of characteristic random processes depending on correlation distance
(b=0.01, 0.1, 1.0, 10.0)
Fig. 3은 추계장의 표준편차 
0.1인 상태에 대하여 상관관계거리별 랜덤 프로세스를 나타낸 것으로서, 상관관계거리에 따른 랜덤 프로세스의 특성을 보여주고 있다. 상관관계거리가 커짐에 따라 랜덤프로세스는 해석영역에서 일정한 값을 가지게 됨을 알 수 있는데, 이론적으로 
인 경우 확률변수상태(random variable state)가 된다.
2.4.1 단일재료-등단면 보에 대한 확정론적 검증해석
기능경사가 없고 등단면인 등질등방성재료 보는 본 연구에서 정식화한 기능경사재료 보에서 식 (1), (2)의 재료 및 기하형상에 대한 두 인수 
를 모두 0으로 택하여 모사할 수 있다. 이러한 가정을 할 경우 식 (3)과 (7)로부터 휨강성은 
로 되며 일반적인 등단면 단일재료 보와 같게 됨을 알 수 있다. 이 경우에 대한 해석을 통하여 확정론적 해석 결과에 의한 하중점 연직변위 및 지지점의 회전변위 모두 보이론과 동일한 결과를 보여, 기능경사재료 보에 대해 기술된 이론 및 정식화의 타당성을 확인하였다.
2.4.2 기능경사재료 보의 응답변화도
보의 재료 및 기하학적 제원으로 탄성계수는 
7.0GPa, 
10.0, 보의 폭은 
1.0m, 
=2.0로 가정하였다. 추계장의 분산계수는 0.1로 가정하였고, 상관관계거리 의 변화에 따른 응답의 분산계수(COV)를 산정하였다. 집중하중이 작용하는 중앙점(x=10.0m)의 연직변위와 
16.0m인 점에서의 회전변위에 대한 분산계수를 제시한다.
대상 기능경사재료 보의 확정론적 해석에 의한 변형형상은 작용하중별로 Fig. 4와 같다.
Fig. 5에서 Y는 연직변위, R은 회전변위를 나타내며, P는 집중하중만 작용한 경우를 P&M은 모멘트하중이 함께 작용한 경우를 나타낸다.
그림에서 보는 바와 같이 제안된 정식화에 의한 결과는 몬테카를로 해석에 의한 결과와 매우 유사한 값을 보여주고 있다. 응답의 분산계수는 중간 상관관계거리에서는 변동을 보이지만, 상관관계거리가 큰 경우 추계장의 분산계수로 접근하는 것을 볼 수 있다. 한 가지 주목해야 할 사항은 응답분산계수가 작용하중의 양태에 따라 동일 자유도에 대하여도 다르게 나타난다는 것이며, 그 값도 추계장의 분산계수보다 클 수 있다는 것이다.
2.4.3 형상인수 
, 재료지수 
의 영향
Fig. 6에서 보는 바와 같이 구조응답의 분산계수는 기하학적 특성에 따라 다른 값을 나타낸다. 보의 좌, 우단에서의 폭의 비() 값이 클수록 응답의 분산계수는 커지지만, 그 증가율이 감소하며 비선형적으로 증가됨을 볼 수 있다. 이는 보의 폭이 커짐에 따라 보의 위치별 강성의 차이가 발생하여 불확실인수의 영향이 증가되기 때문이다.
그러나 식 (1)의 재료지수 에 따르는 두께방향으로의 탄성계수 변화비율인 중립면과 최외곽면의 탄성계수의 비 
의 값은 응답 분산계수에 영향을 주지 않는 것으로 나타났다. 즉, 
이 변화되어도 변위의 분산계수는 동일한 값을 나타낸다. 
의 변화에 따라 구조응답의 1, 2차 모멘트는 크게 변화하지만, 그 비율로 결정되는 분산계수는 동일한 값을 가지게 된다.
2.4.4 추계장 분산계수 
의 영향
추계장 분산계수가 응답변화도 분산계수에 미치는 영향을 보기 위하여 추계장분산계수를 변화시키며 P&M의 하중조건하에서 수행하였다.
Fig. 7에서 보는 바와 같이 본 연구에서 제시한 해석법은 구조응답의 분산계수가 추계장 분산의 값에 따라 선형증가하지만, MCS의 경우 비선형적으로 증가함을 알 수 있다. 이는 식 (10)으로 표현된 탄성계수가 제안 해석법에서는 선형적으로 고려되는 반면, MCS의 경우 직접적으로 역수로 개입되어 
로 구현되기 때문이다. Noh and Park(2006)는 이러한 현상을 모사할 수 있는 등가추계장에 대한 연구를 제시하고 있다. 여기서, 추계장 분산계수를 0.25로 제안한 것은 그 이상의 분산계수에는 랜덤 프로세스의 값이 
-1.0인 경우가 발생되며, 따라서 식 (10)의 가정에서 탄성계수가 물리적인 의미를 가지지 않게 되기 때문이다.
2.4.5 제안 해석법의 효율성
전술한 바와 같이 몬테카를로 해석은 기 설정된 추계장에 대한 통계적 특성치(예, 평균, 표준편차, 추계장의 변동 진동수 등)를 만족하는 랜덤 프로세스(1차원) 또는 랜덤 필드(2, 3차원) 샘플을 수치적으로 생성하고, 이에 대한 확정론적 해석을 수행한 후, 그 결과들에 대한 통계처리 과정으로 응답변화도를 얻는다.
본 연구에서 사용한 스펙트럼 모사법의 경우 랜덤 프로세스의 생성에 소요되는 시간은 생성샘플의 개수, 수치생성에서 사용하는 cosine 함수의 개수인 식 (22)의 N, 그리고 상관관계거리 
가 짧을 경우 aliasing 효과를 제거하기 위한 고려사항 등이 포함된다.
제안 해석법은 하나의 상관관계거리에 대한 해석 시 0.52초의 해석시간이 소요되었다. MCS의 경우 10,000개의 샘플에 대한 해석시간은 8.11초가 요구되었으며, 따라서 제안 방법과 동일한 해석시간에 약 641개의 샘플에 대한 해석을 수행한 것에 해당된다. 일반적으로 MCS에서 적절한 결과를 얻기 위해서는 수천 개의 샘플에 대한 해석이 요구되며, 따라서 해석시간상에서 제안 해석법이 매우 효율적임을 알 수 있다.
샘플의 생성에 소요되는 시간은 10,000개의 샘플 생성시, N=512에서, 34.36초가 소요되어, 이를 함께 고려한다면 제안 해석법의 효율성이 MCS 대비 매우 높음을 알 수 있다.
3. 결 론
본 연구에서는 변단면 기능경사재료 보에 대한 추계론적 유한요소해석을 위한 정식화를 제시하였다. 구조의 불확실성은 보의 축을 따라는 중앙면에서의 탄성계수에 적용하였다. 따라서 중앙면 탄성계수는 축방향으로 공간적 불확실성을 가지게 되며, 두께 방향으로는 확정적 재료지수에 따라 탄성계수가 분포한다.
기능경사재료 보의 추계론적 유한요소해석을 위하여 Hermite 형상함수에 기초한 2-절점, 4-자유도 보 유한요소를 유도하여 사용하였다. 이 경우 개별 보 요소에는 추계적분으로 정의되는 확률변수가 6개 존재하게 되며, 이 확률변수의 평균을 중심으로 하여 변위에 대한 선형 Taylor전개를 사용하여 응답변화도 산정을 위한 1, 2차 모멘트식을 제시하고 이를 이용하여 응답의 분산계수를 산정, 제시하였다.
제안 정식화의 검증을 위하여 스펙트럼모사법을 적용한 MCS 해석을 수행하여 그 결과를 비교하였다. 산정된 응답변화도는 제안한 해석법과 MCS에서 매우 유사한 값을 나타내어 제안 해석법의 타당성을 알 수 있었다. 보의 형상인수, 추계장의 분산계수 등이 구조의 응답변화도에 미치는 영향을 고찰하였고, 해석에 소요되는 시간을 비교한 결과 제안한 해석법은 MCS 대비 그 정확성과 효율성이 매우 높은 것으로 나타났다.
본 연구는 추후 두께방향으로의 탄성계수의 불확실성을 고려하는 연구 및 기능경사재료 평판이나 쉘구조에 대한 추계론적 해석연구의 기초가 될 수 있을 것으로 기대된다.
