2.1 Mullins effect의 적용

고무 혹은 폴리우레탄 등의 초탄성(hyper-elastic) 또는 점탄성(visco-elastic) 재질의 응력을 구하는 수치해석에서는 손상변수(damage variable) $\eta$를 이용하여 변형률에너지 밀도 함수(strain-energy density function) $W(\lambda_{i})$를 변형하여 사용한다.

단축 인장 또는 압축일 경우 Cauchy stress $\sigma_{1}$은

이므로 손상변수(damage variable) $\eta$를 고려할 때, Mullins 모델을 다음과 같이 모델링 한다.

여기서, $\tanh(x)$는 단순한 curve fitting의 역할만을 할 뿐이다. 고무나 폴리우레탄 같은 소재의 강성은 일반적으로 비선형성을 가지고 있고, 이러한 소재들은 대부분 변형이 크다. 일반적으로 대변형 해석에 있어서는 선형계인 Hook’s law를 적용할 수 없으며, Helmholtz free energy $\rho\Psi$에서 변형률에너지(strain energy)만을 추출하여 표현하게 된다.

제 2차 Piola-Kirchhoff 응력(the 2nd Piola-Kirchhoff stress) $\sigma^{PK2}$는 다음과 같이 나타내진다.

여기서 $E_{i}$는 Green의 주변형률(principal strain)이며, $W=W$$(\lambda_{1},\:\lambda_{2},\:\lambda_{3})$이다. 또한, $\lambda$는 연신율(stretch ratio)이라 하고, $\lambda =1+\epsilon_{eng}=1+\dfrac{\Delta L}{L}$이다. 연쇄법칙(chain rule)에 의하여 $\sigma^{PK2}$는 다음과 같이 된다.

한편, Green의 주변형률과 연신율과의 관계는 다음과 같다.

따라서, $\dfrac{\partial\lambda_{i}}{\partial E_{i}}=\dfrac{1}{\lambda_{i}}$이고, 즉, $\sigma_{i}^{PK2}=\dfrac{1}{\lambda_{i}}\dfrac{\partial W}{\partial\lambda_{i}}$가 된다.

Principal Cauchy stress $\sigma_{i}=\lambda_{i}^{2}\sigma_{1}^{PK2}$이므로 $\sigma_{i}$는 다음과 같이 된다.

변형률에너지 함수 $W(\lambda_{1},\:\lambda_{2},\:\lambda_{3})$ 또는 $W(I_{1},\:I_{2},\:I_{3})$의 일반화된 형태는 다음과 같다.

여기서, Invariants $I_{1},\: I_{2},\: I_{3}$ 는 식 (11)~ 식 (13)과 같다.

따라서,

이 되며, 식 (14)에서의 각각의 미분 항들은 다음과 같다.

여기서, $D=\dfrac{1}{2}E_{b}$이며, $E_{b}$는 체적탄성계수(bulk modulus)이다.

식 (14)로부터 가장 일반적으로 사용되는 모델은 $C_{10},\:$ $C_{01},\:$$C_{20}$를 사용하는 Signiorini 모델이다. 이 모델을 사용하면 Cauchy stress $\sigma_{i}$는 다음과 같다.

식 (19)~ 식 (21)에서 $D$와 체적변형(volume change) $\epsilon_{vol}$은 volumetric test로 구할 수 있는데, 본 논문에서는 이 항을 제거하기 위한 방법으로 $\sigma_{2}혹은\sigma_{3}=0$를 적용하였다. 따라서

이 되고, Engineering stress $\sigma_{i,\:eng}=\lambda_{i}\sigma_{i}$ 이므로,

이 된다. 단축시험인 경우에는 $\lambda_{1}=\lambda ,\:\lambda_{2}=\lambda_{3}=\dfrac{1}{\sqrt{\lambda}}$이므로

이 된다. 한편, 강성 $E$는 다음과 같이 계산될 수 있다.

초기위치($\epsilon_{eng}=0$)에서 $\lambda =1$이므로, $E_{0}=6(C_{10}+C_{01})$, 초기전단탄성계수(shear modulus) $G_{0}$는 $G_{0}=2(C_{10}+C_{01})$가 된다.

2.2 응력이완(stress relaxation)

일반적으로 소재의 응력 이완 실험에 있어, 변형률 변화량(strain rate)이 증가하면 구조 감쇠의 영향에 의하여 강성이 증가하게 되는데, 이를 그 소재의 동강성(dynamic stiffness)이라고 한다. 고무나 폴리우레탄 같은 폴리머(polymer) 재료의 동강성과 점탄성 성질을 수학적으로 모델화하기 위해, 다양한 형태의 모델을 이용하게 된다. 본 논문에서는 Prony-series모델(spring-dash pot이 직렬 연결되고, 스프링 하나가 다시 병렬 연결되어 무한히 계속되는 모델)을 이용하였고, Prony-series 모델은 일반화된 Maxwell 모델로 간주하는데, 그림 1에 나타난 바와 같다.

Prony는 시간의 함수로 나타나는 소재의 동강성을 다음과 같이 표현하였다.

여기서, $E_{\infty}$는 $t\mapsto\infty$일 때의 강성을 의미한다. $E_{R}(0)$를 시간이 영(zero)일 때의 강성을 나타낸다면, $E_{R}(t)$는 다음과 같이 다시 표현될 수 있고, 초기 강성을 알고 있을 경우는 더욱 편리하다.

식 (28)을 식 (29)와 같이 다시 표현하고, 정규화된 동강성 $e_{R}(t)$로 표현하면 식 (30)과 같다.

프로니 급수(Prony-series)는 시간과 더불어 감소하다가 일정한 값으로 수렴하는 쇠퇴하는 거동을 나타내며, 이를 수치적으로 근사화 시키기 위해 사용되는 급수이다. 이는 크기가 순차적으로 변하는 시간함수 또는 시간 상수를 가진 유한개의 지수함수들의 조합으로 표현되는 시간에 대한 급수 함수로 정의된다. 또한, 각각의 지수 함수에 곱해지는 상수 $e_{i}$와 $\tau_{i}$는 점탄성(viscoelasticity) 성질을 지닌 재료의 실제 응력이완 현상을 수치적으로 근사화 시키기 위해 효과적으로 사용되고 있다. 한편, 상수 $e_{i}$와 $\tau_{i}$는 응력이완 실험값으로부터, 최소자승법을 적용하여 구할 수 있다.

변형률과 강성이 시간에 따라 변하는 경우의 응력은 다음과 같이 표현된다.

$\dot\epsilon(t)=0$가 되는 경우는 $\epsilon =const$가 되고, $\sigma(t)=E_{R}(t)\epsilon$ 이 된다. Prony-series 모델을 이용하여 표현하면 다음과 같다.

식 (32)를 부분적분하면 식 (33) 또는 식 (34)로 나타내진다.

여기서, $\sigma_{0}(t)$는 시간 $t$에서의 순간 응력을 의미한다.

본 논문의 수치해석 과정은 크게 두 과정으로 나눌 수 있다. 첫째, 2장의 이론에서 언급한 Mullins 효과를 포함한 비선형 응력-변형률 모델의 상수 값들 $C_{10},\: C_{01},\: C_{20}$는 모델과 실험데이터를 사용하고, pseudo-inverse matrix를 이용하여 결정한다. 이들 상수로부터 설계탄성계수 $E^{*}$와 스프링 상수 $k$를 구할 수 있다. 둘째, 2장의 이론에서 언급한 응력이완 부분에 있어, 점탄성 및 동강성을 표현할 수 있는 프로니 급수의 상수들 $e_{i}$와 $\tau_{i}$를 구하는 과정으로 $e_{i}$와 $\tau_{i}$는 최소자승법(least square method)을 이용하여 결정하며, 상용 소프트웨어 MATLAB의 함수인 ‘lsqcurvefit’를 이용하였다.

3.1 스프링 상수와 설계 탄성계수 결정

원통형 시편과 같은 경우, 공학적 관점에서 스프링상수 $k$는 일반적으로 힘과 변위와의 관계로부터 다음으로 나타내진다.

여기서, $\delta$는 압축시험으로부터의 압축량, $A$와 $L$은 각각 소재의 단면적과 길이를 나타낸다. 또한, $E^{*}$는 설계탄성계수 또는 유효강성(effective compressive modulus)으로 시험편의 형상계수(shape factor) S와 압축계수(compressive coefficient) $\phi$의 함수로 다음과 같이 표현된다.

여기서, 형상계수 $S=\dfrac{하중을 받는 단면적}{자유 단면적}$으로 본 논문의 시편의 경우, 다음과 같다.

한편, 압축계수 $\phi$는 표 1에 제시되어 있다.

고무나 폴리우레탄 등은 비압축성(incompressible) 재료로 프와송비(Poisson’s ratio) $\nu$는 $\nu\approx 0.5$에 근접한 값을 갖는다. 따라서 영률(Young’s modulus) $E$는 $E=2(1+\nu)G\approx 3G$이며, 체적탄성계수(bulk modulus) $E_{b}$는 다음과 같이 표현된다.

예를 들어 $E_{b}=1000G$인 경우, $\nu =0.4995$이 된다. 이와 같이 $G,\: E,\: E_{b}$의 $\nu$의 값에 의해 결정되나 $\nu$ 값의 미세 변화가 물성치의 큰 오차를 유발하므로 매우 신중하게 처리되어야 한다.

3.2 응력이완 수치해석

Prony series 모델과 정규화된 강성함수(normalized modulus function $e_{R}(t)$)의 상호관계를 정의하기 위해서 그림 2와 같은 모델을 적용한다.

그림 2의 운동방정식은 그림 3과 같은 일반적인 힘 $f$가 작용하는 경우의 모델에 대한 운동방정식을 표현한 후, 표현하기로 한다. 그림 3모델의 운동방정식은 우선 2개의 식 (39)와 식 (40)으로 표현될 수 있다.

식 (39)와 식 (40)에서 $\dot x_{2}$를 소거하면 다음과 같은 운동방정식이 된다.

식 (41)의 운동방정식과 같이 그림 2모델에 대한 운동방정식도 다음과 같이 표현될 수 있다.

응력이완 중에는 식 (42)에서 $\dot\epsilon(t)=0$ 이므로 식 (42)는 다음과 같이 압축된다.

식 (43)의 해는 제차방정식의 해 $\sigma_{h}(t)=\sigma_{0}\exp(- E_{1}/ C_{1}t)$와 특수해 $\sigma_{p}(t)=E_{\infty}\epsilon_{0}$의 합에 의해 다음과 같다.

그림 2의 Prony-series의 단순 모델을 이용할 때, 식 (30)에서 $i=1$인 경우의 정규화된 강성함수 $e_{R}(t)$는 다음과 같이 표현된다.

식 (44)를 식 (45)와 비교하면 $\tau_{1}=C_{1}/E_{1},\:$$e_{1}=\dfrac{\sigma_{0}}{\sigma_{0}+E_{\infty}\epsilon_{0}}$이 된다. 즉, $e_{1}=\dfrac{\sigma(0)-\sigma(\infty)}{\sigma(0)}$를 의미한다.

4.1 실험장치

고무와 폴리우레탄 같은 폴리머 소재의 정적 압축시험과 응력이완 시험을 위해 다음과 같은 실험장치를 구성하였으며, 실험장치의 사진과 개략도는 각각 사진 1과 그림 4에 나타난 바와 같다. 또한, 시편의 형상 및 치수는 ASTM D695 또는 ISO604에서 제시하는 바와 같이 지름 1/2 inch, 높이 1 inch의 원통 형상이다. 시험에 사용된 고무 시편의 예시는 사진 2에 나타나 있다.

4.2 실험방법

정적 압축시험은 재료의 감쇠효과를 최소화하기 위하여 매우 느린 속도로 수행하며 최대 압축량은 30% 및 50%로 실험하였으며, loading-unloading cycle을 10회 수행하였다. 사진 3과 사진 4는 압축시험을 위해 시편을 장치의 치구에 결합한 상태를 보여준다. Relaxation test는 물질의 감쇠특성을 확인하기 위한 시험으로 시편을 압축 중에 정해진 strain 상태에서 변위를 유지하면 compressive stress가 시간이 지남에 따라 감소하게 된다. 이 감소 상태로 감쇠특성을 정량화한다. 시험은 $\epsilon_{T}$ = 10%, 20%, 30%, 40% 및 50%의 5단계로 수행하였으며 압축 중 1차로 $\epsilon_{T}$ = 10%가 되면 압축을 멈추고 10초간 유지한다. 10초간 stress가 감소하는 상태를 측정한 후 다시 압축을 진행하여 $\epsilon_{T}$ = 20%가 되면 멈추고 10초간 유지한다. 이 과정을 $\epsilon_{T}$ = 50%가 될 때까지 수행한다.

5.1 정적 압축 시험결과

그림 5는 진변형률(true strain)이 최대 30%일 때, 다양한 경도(duro 40, 50, 60, 70)를 갖는 고무 소재의 정적 압축시험을 통해 얻은 공칭응력-변형률 선도이다. 이 경우, 10cyclic test를 행하였다. 이 그림에서 알 수 있는 것은 duro 경도가 높을수록 공칭응력-변형률 선도의 기울기가 큼을 알 수 있으며, 최대 진변형률이 30%일 경우, 식 (1)을 이용한 계산과 같이 공칭 변형률은 약 0.2592 정도임을 알 수 있다.

그림 6은 그림 5의 다양한 경도를 갖는 고무 소재 중, duro 70의 고무 소재에 대해 최대 진변형률 30%일 때 3 cyclic test로부터 얻은 공칭응력-변형률 선도이다. 이 경우, 첫 번째 cycle test(검은색 실선)후, 변형률은 원래의 시작 위치로 돌아가지 못하고 있음을 알 수 있는데 이는 이미 이론에서 기술한 바와 같이 ‘Mullins effect’를 보여주고 있는 것이다. 그 후, 2~3 cyclic test에서는 거의 반복되는 수렴 곡선을 보여주고 있다. 아주 짧은 시간의 충격이 가해지는 시스템에 대한 용도로 본 소재를 사용하는 경우는 설계탄성계수와 스프링 상수를 구하는 데 있어 첫 번째 cycle test 결과의 곡선을 사용할 수 있겠지만, 일반적인 진동문제의 소재로 사용할 경우에는 2 cycle 이후의 수렴되는 곡선을 사용하는 것이 타당할 것으로 사료된다.

그림 7과 8은 최대 진변형률(true strain)이 각각 30%와 50%일 때, 다양한 경도(duro 70, 80, 90)를 갖는 폴리우레탄 소재의 정적 압축시험을 통해 얻은 공칭응력-변형률 선도를 보여준다. 이 경우도 10cyclic test를 행하였으며, 이 그림들 각각의 경우에서도 duro 경도가 높을수록 공칭응력-변형률 선도의 기울기가 큼을 알 수 있으며, 최대 진변형률이 30%일 경우, 식 (1)을 이용한 계산과 같이 공칭 변형률은 약 0.2592, 50%일 경우 공칭변형률은 약 0.3935를 나타내 보이고 있다. 그림 7의 경우는 변형률 변화율이 $\dot\epsilon_{T}=$0.0078/s, 그림 8의 경우는 $\dot\epsilon_{T}=$0.01/s이다.

그림 9는 진변형률(true strain)이 최대 30%일 때, 1st cycle은 배제한 상태의 고무 소재의 공칭응력-변형률 선도이다. 그림 10은 그림 9의 조건과 같지만, 그림 9의 시험한 소재를 1시간이 경과한 후 재시험해서 얻은 공칭응력-변형률 선도이다. 그림 9와 10의 비교에 있어서, 다양한 경도 그래프 모두에 대해 아주 경미한 기울기의 변화는 있지만 큰 변화는 없음을 알 수 있었다. 많은 문헌에 의하면 고분자 소재는 하루 또는 수일이 지나면 1st cycle 후의 소성변형이 회복된다고 나타내고 있지만, 본 논문 실험에서와 같이 1시간 정도에서의 변화는 크게 없는 것으로 나타났다.

그림 11은 duro 90인 폴리우레탄 소재의 3 cyclic test 결과를 보여주고 있다. cycle별로 분리하여 보면 1차 압축 시(loading)에는 높은 강성(검은 선)을 가지며, 복귀 시(unloading)에는 강성이 감소(softening)하는 것을 알 수 있으나, 2차 이후의 압축 시에는 강성이 감소된 상태에서 다시 압축되며(loading) 복귀 경로는 1차와 거의 같은 경로를 따른다는 것을 알 수 있다.

그림 12는 duro 70의 폴리우레탄 소재를 최대 진변형률 30%로 압축 시험한 결과인데, 이 경우 1차 압축 시에만 Mullins 효과에 의한 강성 증대효과가 나타나며, 2번째 이후의 cycle에서는 안정된 상태를 유지한다. 2번째 cycle 이후의 압축 데이터를 사용하여 비선형 모델에 대한 계수 $C_{10},\: C_{01},\: C_{20}$ 구하고, 이 계수는 시험치와 모델을 이용하여 pseudo-inverse matrix를 이용하여 결정한다. 푸른 점선의 기울기는 스프링 상수(engi- neering stiffness) $k= E^{*}A / L$이며, $E^{*}=\Delta\sigma_{eng}/\Delta\epsilon_{eng}$가 되어 시험 결과 값으로 계산할 수 있다.

그림 13과 14는 그림 12의 2번째 압축 곡선을 비선형 모델로 계산하여 얻은 실험값과 수치해 모델의 응력-변형률 선도와 변형률-탄성계수 곡선이다. 그림14의 계수는 $C_{10}=1.3701,\:$ $C_{01}=-0.3223,\:$ $C_{20}=0.0565$이고, 그림 15에서 $k$와 $E^{*}$값은 각각 $k=41.7k N/m,\:$ $E^{*}=8.07MPa$로 나타났다.

그림 15의 고무의 경우도 폴리우레탄 소재와 모두 유사한 경향(cycle 별로 분리하여 보면 1차 압축 시(loading)에는 높은 강성(검은 선)을 가지며, 복귀 시(unloading)에는 강성이 감소(softening)하는 것)을 나타내고 있고, 이후의 그림들에서도 알 수 있듯이, duro 경도에 무관하게 유사한 경향을 나타냄을 알 수 있다. 그림 16은 duro 40의 고무 소재에 있어, 2nd cycle을 기준으로 했을 때의 모델 값과 실험값을 비교한 공칭응력-변형률 선도이고, 그림 17은 이에 따른 수치해 모델의 탄성계수-공칭변형률 그래프를 보여준다. 그림 18은 duro 40의 고무 소재에 있어, 1st cycle을 기준으로 했을 때의 모델 값과 실험값을 비교한 공칭응력-변형률 선도이고, 그림 19는 이에 따른 수치해 모델의 탄성계수-공칭변형률 그래프를 보여준다. 그림 16에서 실험값과 수치해 모델 값의 경미한 차이는 실험에서는 1차 압축 후, 복귀하는 과정에서 일시적인 소성변형이 있었기 때문이며, 그림 18과 같이 1st cycle을 기준으로 했을 때는 수치해 모델과 실험값이 잘 일치하고 있음을 알 수 있다. 실제 설계에는 2cycle 기준에 맞추는 것이 타당할 것으로 사료된다. 그림 16에서 상수 값은 $C_{10}= -0.0642 ,\:$ $C_{01}= 0.3412 ,\:$ $C_{20}= -0.1817$인 값으로 나타났다. 한편, 그림 18에서는 $C_{10}= 0.8273 ,\:$ $C_{01}= 0.4396 ,\:$ $C_{20}= 0.1647$이다.

5.2 응력이완 결과

고무 및 폴리우레탄 소재의 점탄성 성질과 동강성 효과를 관찰하기 위해 응력이완 시험을 하였으며, 이에 대한 결과들은 다음과 같다. 그림 20은 듀로 경도 70의 폴리우레탄 재료의 relaxation test를 한 결과이다. 검은 선은 strain을 붉은 선은 stress를 나타낸다. 그림에 나타낸 바와 같이 strain을 증가시키다 멈춘 후 일정하게 유지시키면서 stress의 감소를 계측하였다. $\dot e_{R}(t)$를 시험적으로 결정하기 위해서 그림 20과 같이 5개 구간에서 순간적으로 $\dot\epsilon(t)=0$가 되는 구간을 만들었다. 각 구간은 $\epsilon_{T}=$10%, 20%, 30%, 40%, 50%가 되는 지점이다. $\epsilon$ 또는 $\epsilon_{T}$가 일정한 구간에서 stress가 감소한다. 이 부분을 나타낸 그림이 그림 21이다. 압축 변형량이 작을수록 감쇠 폭이 커짐 ($\tau_{i}=\eta_{i}/E_{i}$가 작아짐)을 알 수 있다. Creep 해석과 같이 장시간을 필요로 하는 시험의 경우에는 매우 큰 (100 s이상) $\tau_{i}$가 필요하지만 그렇지 않은 충격문제 등에서는 매우 짧은 시간의 반응만 필요하므로 긴 시간은 무시한다. 시험 데이터를 이용하여 $e_{R}(t)$의 매개변수 $e_{i}$ 와 $\tau_{i}$는 최소자승법을 이용하여 결정하였으며, MATLAB의 함수인 lsqcurvefit를 사용하였다.

표 2는 duro 경도 70의 폴리우레탄 소재의 응력이완이 되는 부분만을 따로 나타낸 그림 21의 결과 그래프로부터 식 (45)의$e_{R}(t)=1-e_{1}[1-\exp(-t/\tau_{1})]$$=(1-e_{1})+e_{1}\exp(-t/\tau_{1})$을 나타낸 식들이다. 즉, 식 (45)의 $e_{1}$과 $\tau_{1}$을 구하여 정리한 식이다. 진변형률 10% ~ 50%에 따른 각각의 식이 나타나 있다.

그림 22와 23은 각각 duro 경도 50의 고무 소재에 대한 응력이완 시험곡선과 정규화된 응력곡선을 $\epsilon_{T}=$10%, 20%, 30%, 40%, 50%에 대해 나타낸 그림들이다. 또한, 그림 24와 25는 각각 $\epsilon_{T}=$10%에 대해, 시간이 매우 짧은 구간(약 0~0.32 s)과 약간 긴 구간(약 0.3~9 s)에 대해 정규화된 응력곡선을 curve fitting한 그래프이다. 표 3은 $\epsilon_{T}=$10% ~ 50%의 각각에 대해 시간이 매우 짧은 구간(0~0.32 s)과 약간 긴 구간(0.3~9 s)에서 duro 50의 고무 소재에 대한 정규화된 modulus함수 $e_{R}(t)$ 식을 나타낸 표이다.