Thermosetting polymers are extensively used in various industries, including shipbuilding, automobiles, and batteries, due to their excellent mechanical properties and chemical resistance derived from their crosslinking network. However, crosslinked networks of polymers are generally considered be irreversibly formed; thus, recently significant attention has been given to covalent adaptable networks (CANs). Still, the mechanisms and dynamics of the bond exchange reaction (BER) are yet unknown, and also the molecular modeling approaches for BER vary across previous studies. Therefore, this study aims to describe BER in crosslinked polymeric systems via nonequilibrium molecular dynamics (NEMD) simulations. Specifically, the stress behavior of these systems is examined by quantitatively analyzing the degree of BER in conjunction with structural and rheological properties. These finding would provide insights into the intrinsic molecular mechanisms and physical origin of CANs.

열경화성 고분자는 가교 구조로 인한 우수한 기계적 특성 및 화학적 저항성 등에 의해 다양한 산업에서 널리 사용되고 있다. 그러나 일반적으로 가교 구조는 비가역적으로 형성되므로, 최근에는 공유 적응형 네트워크가 주목받고 있다. 하지만 결합 교환 반응의 메커니즘과 동역학은 여전히 명확하게 규명되지 않았고, 분자 모델링 접근법 역시 이전 연구마다 다르다. 따라서 본 연구는 비평형 분자동역학 시뮬레이션을 통해 가교 고분자 시스템 내 결합 교환 반응을 모사하고자 하였다. 특히, 구조적 및 유변학적 특성과 관련하여 교환 반응 정도와 함께 정량적으로 분석함으로써, 가교된 고분자 시스템의 응력 거동을 규명하였다. 이러한 발견은 CAN의 분자 메커니즘 및 물리적인 기원에 대한 중요한 통찰을 제공할 수 있을 것이다.

This study reveals that probability of bond exchange reaction, along with the frequency and lifetime of bond creation and breaking, significantly influence the stress relaxational behaviors in crosslinked polymers, via nonequilibrium molecular dynamics (NEMD).

Keywords: molecular dynamics simulation, stress relaxation, covalent adaptable networks, dynamic covalent bonds.

이 성과는 정부(과학기술정보통신부)의재원으로 정보통신기획평가원의 지원(2020-0-01791)과 한국연구재단의 지원(RS-2022-01165711)을 받아 수행된 연구임. 또한 본 연구는 2024학년도 동의대학교 교내연구비에 의해 연구되었음(202401330001).

저자는 이해상충이 없음을 선언합니다.

열경화성 수지는 가교 네트워크 특성으로 인해 기계적 강도 및 열 안정성, 화학적 저항성 등이 우수하다고 알려져 있다.¹ 일반적으로 고분자 가교로 인한 네트워크 형성은 비가역적이기 때문에, 재가공이나 재활용이 어렵다. 그러나 최근 가역적 네트워크를 형성할 수 있는 에폭시 고분자 및 가교제가 개발되면서, 친환경적이며 재활용 가능한 가교 고분자 소재에 대한 관심이 커지고 있다. 가역적 네트워크 형성을 가능하게 하는 동적 공유 결합(dynamic covalent bonds, DCBs)은^2,3 전통적인 공유 결합과 달리 온도, UV 등의 특정 조건 아래에서 상대적으로 낮은 활성화에너지를 가지면서 가역적으로 결합 생성, 교환 및 끊어짐이 가능하며, 이를 통해 형성된 시스템을 공유 적응형 네트워크(covalent adaptable networks, CANs)라고 한다(Figure 1). 따라서 기존의 열경화성 수지의 특성을 유지하면서 자가 치유, 응력 완화, 재활용성 등의 흥미로운 거동을 보인다. 예를 들면, Scott 등에⁴ 따르면 빛의 조사에 따라 응력 완화를 촉진시키는 결합 교환 반응(bond exchange reaction, BER)이 가능함을 확인하였으며, Leibler 등에⁵ 의해 고분자 reptation 모델의 sticker를 도입하여, 동적 네트워크의 응력 완화 모델 이론을 해석하는 등의 결과들이 보고되었다.
특히 Leibler 등에⁵ 의해 발견된 비트리머(vitrimer)는 associative CAN의 하나로, 유리전이온도(T_g) 외에 topology freezing transition temperature(T_v)라고 알려진 다른 전이 온도 특성을 보인다. 비트리머는 T_v의 온도에 도달하면 교환반응이 활발해지면서 점도 거동이 바뀌기 때문에 점탄성 유체의 성질을 보이게 된다. 그렇기에 T_v에 영향을 주는 인자 및 메커니즘에 대한 심도 있는 이해가 필요하나 명확한 T_v의 정의는 아직 정해지지 않았기 때문에, 다양한 실험적 접근 방식을 기반으로 측정되고 있다.^6-8 현재까지 많이 연구된 DCBs가 가능한 소재들은 disulfide bond, imine bond, ester bond, Diels-Alder (DA) bonds 등을 포함하고 있다고 알려져 있으며,⁹ 이외에도 다수의 관련 연구가³ 진행되어왔지만 교환반응에 의해 네트워크가 재배열되고, 부과된 변형이 완화될 때의 미시적인 메커니즘에 대한 면밀한 해석은 아직 부족하다.
분자 동역학(molecular dynamics, MD)은 뉴턴역학에 기반하여 분자의 운동을 서술하는 계산 기법으로, 실험에서 관찰하기 어려운 미시적인 분자 규모의 역학 및 물성을 탐구할 수 있다. 따라서 본 기법을 바탕으로 결합 교환 반응을 모사하여 CANs의 메커니즘 및 그에 의한 물성을 조사할 수 있다. 예를 들어, Hua 등은¹⁰ PCFF forcefield를 사용한 분자동역학 계산을 통해 에폭시 수지 DGEBA로 형성된 가교 네트워크의 거시적인 구조 뿐만 아니라, 가교 부위 간 거리, 사슬 각도 등의 미시적인 속성을 조사하여 scaling 이론과 비교 분석하였다. 또한 Singh 등은¹¹ 교환 반응이 일어날 확률을 사용하여, T_v 근처의 온도에서 uniaxial tension에 의한 stress가 완화되는 것을 CVFF forcefield로 원자단위 모델링을 통해 모사하였다. 특히 교환 반응으로 인한 자가 치유 특성을 직접 관찰하였고, 온도에 따라 mobility 증가로 인한 교환 반응의 활성화를 확인하였다. 추가적으로 Rottach 등은¹² FENE 모델을 활용하여 가교된 고분자 네트워크의 응력 완화 거동을 계산하였고, slip tube model로 정확하게 예측됨을 확인하였다. 이외에도 free surface를 두고 교환 반응에 의해 표면 welding을¹³ 예측하거나, accelerated ReaxFF 알고리즘을 활용하여¹⁴ uniaxial tension에 의한 탄소섬유 강화 고분자 복합재의 기계적 물성을 계산하는 등의 선행 연구들이 있다. 하지만 여전히 교환 반응의 분자 메커니즘에 대한 이해는 아직 부족하며 특히 교환 반응에 의한 네트워크 재배열과 응력 완화 거동 간의 관계성 및 미시적인 정보는 부족한 상태이다.
본 연구에서는 에폭시 단량체인 polyethylene glycol diglycidyl ether(PEGDGE)와 crosslinker인 2-aminophenyl disulfide(AFD)가 가교된 시스템에 uniaxial tensile을 가한 후, 응력 완화 거동을 비평형 분자동역학 시뮬레이션을 통해 조사하였다. 또한 효과적인 동적 교환 결합을 모사할 수 있는 전산모사 방법론을 다루고 DCBs를 온도에 대한 확률 모델로 묘사하여 결합 생성, 교환 및 파괴를 구현하였으며, 교환 반응의 메커니즘을 구체화하여 그에 따라 응력 완화 거동이 어떻게 변화하는지를 분자레벨에서 연관 지어 분석하였다.

Modeling of Crosslinked System. 에폭시 monomer PEGDGE (296 mol, n=15, 792 g/mol)와 crosslinker 2-AFD(148 mol, 248 g/mol)를 2:1로 섞은 cell을 구축하고, BIOVIA사의 Materials Studio 2024를 활용하여 300 K, 1atm에서 NPT ensemble로 충분히 섞일 수 있게 계산을 수행한 후, NVT ensemble을 통해 구조를 안정화하여(73.31 × 73.31 × 73.31 Å) dimension의 1.14 g/cm³ density를 갖는 시스템을 얻었다(Figure 2).
이 후 Figure 2에 표시된 R1과 R2 간 crosslinking을 만들기위해, LAMMPS를¹⁵ 사용하여 crosslinking density가 0.9에 도달할 때까지 NPT ensemble아래 시뮬레이션을 수행하여 crosslinked PEGDGE 시스템을 얻었다. 또한 가교가 완료된 시스템에 대해서 NPT ensemble로 충분한 시간동안 relax를 시켜주었다. 각 분자들의 상호작용은 PCFF forcefield로¹⁶ 계산되었으며, 비 결합 상호 작용은 분산 및 정전기력에 대해 각각 2.5 s 및 Ewald summation의¹⁷ cut-off 거리를 갖는 Lennard-Jones(LJ) 12-6를 사용하여 평가되었다.



여기서, r은 원자 간 거리, σ와 ε은 각각 구조 및 에너지 parameter를 의미한다.
Molecular Dynamics Simulation with BER. Crosslinking이 완료된 PEGDGE에 대해 LAMMPS에 기반하여, 300 K 아래 2×10⁹/s 속도로 uniaxial tension을 가하고 strain이 0.3에 도달하였을 때 tension을 제거하는 응력 완화 NEMD 시뮬레이션을 수행하였다. 응력이 완화되면서 발생하는 disulfide bond의 결합과 해리를 온도에 따른 확률함수로 서술하고,¹¹ BER 발생 후 시스템을 이완시키기 위한 MD cycle(NPT-NVT ensembles)을 같이 진행하였다.¹⁸ 각각의 MC 단계에 대해, 기존의 가교 결합은 1 Å의 cuf-off 거리 이상일 때, 9 Å 거리 이하일 때 온도 의존적 확률에 따라 가교 결합이 제거 및 생성된다. 결합 생성과 해리 단계는 20 ps씩 진행되었으며, 각 단계 후에는 시스템 내 결합 정보 변경에 따른 불안정한 에너지를 완화시키기 위해 50 ps 이상의 MD cycle을 진행하였다. 여기서 사용된 단계별 시간은 BER 발생과 energy minimization이 충분하게 이루어질 수 있는 시간으로 설정하였고, 모든 단계에 대해 100번의 iteration을 진행하였다. 특히 온도에 따른 부피 변화를 계산하여 T_g를 산출하고 그에 맞는 T_v 온도를 가정하여 확률 모델을 도출하였다. 시뮬레이션의 timestep은 1 fs로 수행되었으며, Nosé-Hoover thermostat^19,20 아래 periodic boundary condition이 적용되었다.

Crosslinked System and BER. Crosslinking이 완료된 PEGDGE의 밀도는 density는 300 K에서 약 1.1 g/cm³로 계산되었으며 uniaxial tension을 가하기 전, 온도에 따른 specific volume을 계산하여 나타내었다(Figure 3). 시스템이 충분히 평형화된 후, T_g를 계산하기 위해 annealing한 후 quenching하는 기법을 사용하였다. 먼저 온도는 500 K까지 올리고 NPT ensemble로 대기압 하에서 500 ps 이상 평형화하였고, 이 후 20 K씩 감소시키면서 100 K까지 계산을 수행하였다. Cooling 시작 전 구조를 충분히 평형화시켰으며, 온도 감소 폭이 20 K이기 때문에 구조 최적화는 빠르며 오차는 작았다. 또한 각 온도의 specific volume 값은 시뮬레이션 계산시간 동안 샘플링한 후 평균하여 도출하였다. Figure 3에 표시한 것과 같이, 특정 온도 근처에서 기울기의 불연속점이 나타나게 되는데 이는 T_g를 의미한다.²¹ 여기서의 T_g는 243.1±5.2 K로, 실험값(223–259 K)의^22-24 범위와 비슷하였다. 그러나 MD 시뮬레이션은 실험과 대비하여 냉각속도가 빠르고, 실험의 시료보다는 crosslinking 간의 간격이 짧기 때문에 예측 값이 높은 경향이 있다. 여기서 계산한 T_g를 바탕으로 T_v의 온도를 약 320 K으로 가정하여, 이에 맞는 BER 확률함수를 식 (2)로 나타내었다.



w와 α는 각각 transition width 및 확률 관련 변수를 의미한다.¹¹ 이 후 평형화 된 crosslinked PEGDGE에 uniaxial tension을 가하면서 strain(γ)에 따른 stress(s) 값을 시스템의 snapshot과 함께 나타내었다(Figure 4).
Strain이 약 0.3에 도달할 때까지 시뮬레이션을 진행하였으며, strain을 가하는 동안은 BER 알고리즘을 적용하지 않았다. 다른 선행 계산 연구와 같이,²⁵ strain이 증가함에 따라 stress가 증가하는 경향이 비슷한 것을 알 수 있다.
Stress Relaxation. strain이 가해진 시스템에 대해서 300, 350, 400 K에 맞는 BER 확률을 계산한 뒤, 3가지 시스템으로 나누어 BER 발생 정도를 달리하였을 때 응력이 어떻게 완화되는지 계산하였다. 먼저 앞서 방법론에서 언급하였던 BER 구현 알고리즘 및 함수를 사용하였을 때, Figure 5와 같이 가교 결합의 생성 및 파괴가 일어나는 것을 확인하였고, 이를 통해 BER 알고리즘이 응력 완화 시뮬레이션 중 잘 적용됨을 알 수 있었다.
Figure 6은 시간에 따라 3가지 시스템의 응력이 완화되는 거동과 system의 configuration을 보여준다. All-atom 계산 특성상 atomistic stress의 fluctuation이 심하므로, block averaging을 거친 후 Figure 6(a)에 나타낸 식 Kohlrausch–Williams–Watts (KWW) 함수로 fitting하여 도출하였다. Figure 6(a)에서 확인할 수 있듯이 3개의 시스템 모두 시간에 따라 응력이 완화되며, crosslinking이 있기 때문에 응력은 완벽하게 완화되지는 않는다. 특히 온도가 증가함에 따라 BER의 확률이 높아지면서, 가역적인 결합의 교환반응으로 인해 변형된 분자 구조들의 재배열이 빨라지고 그에 의한 완화가 촉진됨을 확인할 수 있다. 또한 KWW 함수로 얻은 relaxation time(τ)을 계산하여 온도 및 BER 확률에 따라 비교해보았을 때, 300 K보다 400 K 시스템의 relaxation time이 약 2배이상 짧았다. 또한 응력 완화 과정 중 특정 시간의 configuration을 살펴보면 실제 시뮬레이션 time은 같음에도 불구하고, 400 K의 경우가 분자 내부 구조들이 좀 더 relax가 되었음을 알 수 있다(Figure 6(b)). 물론 온도가 증가함에 따라 분자들의 thermal Brownian motion이 증가하면서 응력 완화되는 속도가 촉진되는 효과도 있으나, 그 역학에 의한 촉진보다 더 빠른 경향을 보였으며 이는 BER에 의해 가교 네트워크가 활발히 변하면서 응력의 완화에 도움을 준 것이라고 볼 수 있다. 여기서 중요한 것은 BER의 메커니즘별로 응력 완화에 어떤 영향을 주는가이다. 앞서 설명했지만 BER은 기존의 가교 결합을 끊거나, free radical들이 새로운 가교 결합을 형성하는 모든 과정을 의미한다. 여기서 분자의 free volume 측면에서 고려하면, 기존의 가교 결합을 많이 끊어버릴수록 free-end가 많이 형성되어 분자의
mobility를 높일 수 있고, 새로운 가교 결합을 만드는 것은 mobility를 낮게 할 수 있다. 또한 lifetime 측면에서 고려하면, 기존의 가교 결합을 자주 끊을수록 그리고 새로이 형성된 가교 결합의 유지 시간이 짧을수록 mobility가 높아질 수 있으므로 응력 완화 촉진에 도움을 줄 것이다. 즉 단순히 BER 반응 자체의 빈도도 중요하지만 결합 생성 및 파괴의 비율과 lifetime 역시 응력 완화에 중요한 변수 중 하나가 될 수 있을 것이다. Figure 7은 각 시스템에 대하여 시간에 따라 발생한 모든 BER 반응의 수(N_BER)를 전체 시스템 내 가능한 crosslinking 수로 정규화하여 나타낸 것이다. 300 K의 경우, 수학적으로 BER 확률이 0은 아니지만 매우 작은 값이기 때문에 BER이 발생하지는 않았다. 그에 반해 400, 450 K의 경우 시간에 따라 누적되는 BER의 수가 증가하는 것을 알 수 있다. 여기서 BER은 결합의 생성 및 파괴를 모두 포함하는 수치이기 때문에, BER의 미세한 메커니즘에 따른 응력 영향을 면밀하게 분석하기에는 정보가 부족하다. 따라서 위의 시뮬레이션은 결합의 생성 및 파괴를 같은 비율로 발생하게 하였으나, 임의로 그 비율을 조정함에 따라 그리고 결합이 생성되었을 때 유지되는 시간 및 결합이 파괴되었을 때 free radical로 유지하는 시간인 lifetime(t_life) 값을 조정하여 기존에 했던 350 K의 시뮬레이션과 비교하여 응력 완화 거동이 어떻게 달라지는 지를 Figure 8에 나타내었다. Figure 8에서 검은색 데이터는 350 K일 때의 응력완화 거동이며 BER의 결합 생성 및 파괴의 비율과 각 t_life을 같게 둔 것으로 Figure 6(a)의 데이터이다. 그리고 파란색의 점선들은 새로이 생성된 가교 결합이 유효하는 시간인 t_life이 결합이 파괴되어 free radical로 존재하는 t_life과 대비한 비율에 따라 계산한 것이다. 그리고 초록색의 그래프는 t_life는 같은 비율로 유지하되, 전체 BER 횟수 중 결합이 파괴되는 횟수의 비율을 늘리면서 계산한 것을 의미한다. Figure 8에서 확연히 알 수 있듯이, BER 자체도 중요하지만 BER 단계 내 lifetime 역시 응력 완화 거동에 중요한 영향을 미칠 수 있으며, 전체 BER 횟수 중 결합의 생성과 파괴의 절대적인 수치도 응력 완화에 주요 인자 중 하나임을 확인할 수 있다. 특히 결합이 파괴되면 생성되는 free radical들이 chain의 mobility를 높여주어, 응력 완화 촉진의 주 원인 중 하나가 될 수 있다. 이러한 요인을 정량적으로 분석하기 위해, Figure 8에서 얻은 relaxation time을 350 K일 때를 기준으로 표준화하였을 때, 각 시스템의 relaxation time을 N_b/N_c와 t_life에 대한 함수로서 표현하였다(Figure 9).
x축명에 표시된 것처럼 두 변수에 대해 power law로 묘사하였고, A와 B는 각각 약 1.5 및 7.5로 산출되었다. 계산적으로 bond breakage와 creation의 비율 및 가교 결합의 평균 유지 step을 조정하여 얻은 결과지만, 각 인자들에 의해 stress relaxation 거동이 영향을 받을 수 있다. 주의할 것은 두 변수 모두 비율이기는 하나, 빈도 및 시간의 단위가 다르므로 절대적인 값 및 power law 상수 크기로만 영향도를 비교할 수 없을 것이다. 다만, N_b/N_c에 좀 더 영향을 받을 것으로 예상할 수 있다. 이는 bond creation의 lifetime이 줄어들수록 새로운 가교결합이 오래 유지되지 않고 free radical로 바뀌며 그에 의한 mobility 상승 및 시스템의 relaxation을 촉진시킬 수 있으나, 그만큼 bond creation할 수 있는 free radical이 많아질 수 있는 약간의 trade-off 경향을 보일 수도 있기에 N_b/N_c에 좀 더 의존할 것으로 판단된다. 마찬가지로 N_b/N_c 역시 lifetime과 비슷한 개념으로 이해될 수 있으므로, N_b/N_c 증가에 의한 free radical의 증가(bond creation의 lifetime 감소) 효과를 계산적으로 분리하는 것도 하나의 중요한 모델링의 요소가 될 수도 있을 것이다. 따라서 BER 모델링을 할 때 기존 선행 연구에서 주로 사용하였던 변수인 BER 확률 혹은 횟수 외에도 각 메커니즘별 lifetime과 발생 비율도 고려해서 모델링을 하여야 좀 더 근본적인 분자레벨의 메커니즘 및 그에 맞는 CANs의 설계가 가능할 것이다.

본 연구는 동적 공유 결합을 갖는 가교 고분자 네트워크의 응력 완화 거동을 비평형 분자동역학 시뮬레이션을 통해 모사함으로써, 동적 공유 결합의 분자 메커니즘 및 응력 완화에 미치는 영향을 규명하고자 하였다. All-atom 모델을 사용하여, 가교 고분자 구조를 형성하고, BER 확률 함수를 대입하여 uniaxial tension이 가해진 crosslinked PEGDGE의 응력 완화를 BER 발생 정도에 따라 계산하였다. BER 반응을 정량적으로 분석한 결과, BER 반응의 확률뿐만 아니라 새로운 가교 결합의 형성과 기존 가교 결합의 파괴의 각각의 빈도와 각 메커니즘의 lifetime 역시 물성에 큰 영향을 끼칠 수 있다는 것을 확인하였다.
기존의 BER 관련 연구에서는 주로 BER이 발생할 수 있는 외부 자극 요소, BER 정도 및 그에 따른 완화 시간 위주로 보았다면 본 연구를 통해 좀 더 미시적인 관점에서의 면밀한 해석이 필요함을 알 수 있다. 또한 CANs를 형성하는 다양한 monomer 및 crosslinker가 있고, 가교 조건이 많고 그에 따라 복합적인 가교 구조를 만들기 때문에 아직 연구되어야 할 부분은 여전히 많으며 본 연구의 결과가 항상 적용될 수는 없을 것이다. 하지만 이와 같은 분자레벨에서의 해석은 BER의 메커니즘과 그에 따른 물성 변화에 대한 이해를 높일 수 있고, 향후 CANs의 설계와 응용 가능성을 제시하는데 기여를 할 수 있을 것으로 예상된다.