# Effect of shear flow on the transverse thermal conductivity of polymer melts

Kotaro Oda<EMAIL_ADDRESS>Shugo Yasuda<EMAIL_ADDRESS>Graduate School of Information Science, University of Hyogo

###### Abstract

The effect of shear flows on the thermal conductivity of polymer melts is
investigated using a reversed nonequilibrium molecular dynamics (RNEMD)
method. We extended the original RNEMD method to simultaneously produce
spatial gradients of temperature and flow velocity in a single direction. This
method enables the accurate measurement of thermal conductivity in the
transverse direction to shear flow.

Our main finding is that the stress-thermal rule (STR) (i.e., the linear
relationship between anisotropic parts of the stress tensor and the thermal
conductivity tensor) holds for entangled polymer melts even under shear flows
but not for unentangled polymer melts. Furthermore, once entanglements form in
polymer chains, the coefficient of the STR (i.e., the stress-thermal
coefficient) remains independent of the polymer chain length.

These observations align with the theoretical foundation of the STR, which
focuses on energy transmission along the network structure of entangled
polymer chains [Brule, Rheol. Acta 28, 257 (1989)]. However, under driven
shear flows, the stress-thermal coefficient is notably smaller than that
measured in the literature for a quasi-quiescent state without external
forces. Although the mechanism of the STR in shear flows has yet to be fully
elucidated, our study confirmed the validity of the STR in shear flows. This
allows us to use the STR as a constitutive equation for computational thermo-
fluid dynamics of polymer melts, thus offering broad engineering applications.

stress-thermal rule, reversed nonequilibrium molecular dynamics, polymer
melts, entanglement

##  I INTRODUCTION

Computer simulations of polymer melts have been actively conducted in the
chemical industry in recent years to achieve more efficient processing and
product design than traditional experimental approaches. However, accurately
predicting the flow behavior and temperature distribution of polymer melts is
still challenging because of their complex thermo-rheological nature. For
example, in high-speed polymer processing, a polymer melt has a nonuniform
temperature distribution due to significant viscous heating, leading to
complicated thermo-rheological behavior. In fact, in our previous studies[1,
2], a multiscale simulation combining molecular dynamics and computational
fluid dynamics demonstrated that polymer melts can exhibit intriguing
transitional behaviors in terms of their stress-optical rule. This behavior
arises from the competition between flow deformation and thermal agitation
induced by significant viscous heating, both of which influence chain
orientation.

Apart from viscous heating, the influence of the anisotropy of thermal
conductivity, an aspect not previously addressed in our previous studies, is
expected to play a crucial role in temperature distribution during polymer
processing. The anisotropy of the thermal conductivity of polymer melts has
recently been studied by various researchers.

Macroscopically, the thermal conductivity tensor $\bm{\kappa}$ is defined by

$\displaystyle\mathbf{q}=-\bm{\kappa}\cdot{\nabla T},$ (1)

where $\nabla T$ is the temperature gradient and $\mathbf{q}$ is the heat flux
vector. For polymer melts, Brule [3] derived a linear relationship between the
thermal conductivity tensor and the stress tensor, which is called the stress-
thermal rule (STR), based on a network theory for polymer melts assuming that
heat is preferentially transmitted along the backbone of a polymer chain.
Specifically, the STR is expressed as follows:

$\displaystyle\bm{\kappa}-\frac{1}{3}\mathrm{tr}(\bm{\kappa})I=\kappa_{eq}C_{\mathrm{t}}\left(\bm{\tau}-\frac{1}{3}\mathrm{tr}(\bm{\tau})I\right),$
(2)

where $\bm{\tau}$ is the stress tensor, $I$ is the identity matrix,
$C_{\text{t}}$ is the stress-thermal coefficient, and $k_{\text{eq}}$ is the
thermal conductivity in the equilibrium state. Although network theory
considers a network structure consisting solely of polymer segments (it
ignores interactions between polymer chains that are not connected to the
network), experimental and simulation studies [4, 5, 6, 7, 8, 9, 10] have
confirmed that the STR indeed holds for the anisotropic parts of the thermal
conductivity tensor and the stress tensor for polymer melts under elongational
deformation. An STR has also been confirmed for polymer melts under shear
deformation obtained after the cessation of shear flow [11, 12, 13, 14].
However, as far as the authors know, an STR has yet to be confirmed for
polymer melts under shear flows. When an STR is confirmed in a flowing system,
one can use it as a constitutive relation for the computational thermo-fluid
dynamics of polymer melts, thus offering broad engineering applications.

Equation (2) indicates that in the uniform shear flow described as
$\nabla\mathbf{v}=\dot{\gamma}\bm{\mathrm{i}}\bm{\mathrm{k}}$ (where
$\mathbf{i}$ and $\mathbf{k}$ are the unit vectors in the $x$ and $z$
directions, respectively), the thermal conductivity in the transverse
direction to the shear flow, $\kappa_{zz}$, decreases since the normal stress
in the transverse direction to the shear flow, $\tau_{zz}$, decreases for
polymer melts. This behavior is completely different from that for simple
fluids. Davis and Evans derived a thermal conductivity tensor for uniform
shear flows of Lennard-Jones (LJ) fluids in Ref. [15], which posits that the
thermal conductivity transverse to shear flows is written as

$\kappa^{\mathrm{LJ}}_{zz}=\kappa_{0}+c\dot{\gamma}^{2},$ (3)

where $\kappa_{0}$ is the thermal conductivity according to Fourier’s law in
the quiescent state and $c$ is a small positive constant. The above equation
indicates that for simple fluids, $\kappa^{\mathrm{LJ}}_{zz}$ is almost
constant when the shear rate $\dot{\gamma}$ is low but increases as
$\dot{\gamma}$ increases when $\dot{\gamma}$ is sufficiently high. Thus, the
STR could significantly contribute to the distinctive behaviors observed in
the thermo-fluid dynamics of polymer melts compared to those of simple fluids.

In this paper, we investigate the STR (i.e., Eq. (2)) of polymer melts under
uniform shear flows using molecular dynamics (MD) simulation. For simple
fluids, the thermal conductivity under shear flows has been calculated by
various methods. Davis and Evans utilized a Green-Kubo method to calculate the
thermal conductivity in shear flow [15]. Smith employed a wall-driven Couette
system with different wall temperatures and determined the thermal
conductivity in shear flow by measuring the heat flux and temperature gradient
between the walls [16]. However, in the case of polymer melts, the
computational time significantly increases due to the slow dynamics of
entangled polymer chains when attempting to calculate the thermal conductivity
using the Green-Kubo method. In the wall-driven approach, the presence of a
wall can significantly alter the molecular structure, leading to results that
may differ substantially from those of bulk characteristics.

To overcome these issues, we employ the reversed nonequilibrium molecular
dynamics (RNEMD) method [17, 18]. This method involves imposing an artificial
flux in the system by swapping the momentum (or kinetic energy) of molecules
at different locations, which in turn generates a shear (or temperature)
gradient as a physical response within the system. Using the RNEMD method to
calculate the thermal conductivity under shear flows has several advantages.
First, it can simultaneously generate a uniform shear flow and a uniform
temperature gradient within an MD cell with periodic boundary conditions. We
modified the original RNEMD method to achieve the simultaneous generation of
shear flow and a temperature gradient. This modification allows us to measure
the thermal conductivity transverse to the shear flow in the bulk of the
polymer melt. Additionally, the RNEMD method conserves the total momentum and
energy of molecules within the MD cell, even when shears and temperature
gradients are applied. This property leads to accurate computations and
enables us to significantly reduce the amount of simulation data compared to
nonconserved systems such as wall-driven systems.

This paper is organized as follows. In Sec. II, we explain our modified RNEMD
method for calculating the thermal conductivity in the direction transverse to
the shear flow. Section III provides numerical results for both simple fluids
and polymer melts. In Sec. III A, we test the validity of our modified RNEMD
method by comparing our numerical results with Eq. (3) for LJ fluids. In Sec.
III B, the relationships between the thermal conductivity, stress tensor, and
bond-orientation tensor are explored, and the STR is discussed (2). For model
polymer melts, we utilize the Kremer-Grest (KG) model [19]. Finally, we
present concluding remarks in Sec. IV.

##  II METHODOLOGY

We applied the RNEMD method to calculate the thermal conductivity in the
transverse direction under shear flow. As shown in Fig. 1, both the heat flux
and the $x$-momentum (i.e., the momentum in the $x$-direction) flux are
imposed in the $z$ direction in a periodic simulation cell with dimensions
$L_{x}\times L_{y}\times L_{z}$.

Figure 1:  Illustration of our RNEMD method. We impose both a nonphysical
momentum flux and a heat flux by swapping momentum and kinetic energy between
a “hot” slab and a “cool” slab in the simulation box. The molecular system
physically responds to these imposed fluxes by generating both a shear flow
and a temperature gradient.

The simulation cell is divided into slabs along the $z$ axis. The central slab
at $z=L_{z}/2$ is assigned to the ”hot” slab, and the edge slab at $z=0$ and
$z=L_{z}$ is assigned to the ”cool” slab. Additionally, the atoms inside the
hot slab are driven in the $+x$ direction, while those inside the cool slab
are driven in the $-x$ direction. The heat flux is generated by swapping the
kinetic energy of velocities $y$ and $z$,
$\frac{m}{2}({v_{y}}^{2}+{v_{z}}^{2})$ ($m$ is the mass of the particle),
between the hot and cold slabs at each constant time interval $\Delta
t_{\text{k}}$. This process involves selecting two atoms: one with the highest
kinetic energy in the cool slab and the other with the lowest kinetic energy
in the hot slab. Then, the $y$ and $z$ components of the velocities are
swapped within the pair. Note that in our approach, we do not swap the $x$
component of the velocities in the kinetic energy (KE) swapping. This differs
from the original RNEMD method. This modification can prevent interference
with the shear field driven in the system by $x$-momentum swapping.

The $x$-momentum flux is also generated by swapping the $x$-momentum between
the hot and cold slabs at each constant time interval $\Delta t_{m}$, where we
select the atom with the largest $v_{x}$ in the cold slab and the one with the
largest $-v_{x}$ in the hot slab and then swap the $x$ components of their
velocities within the pair.

At the macroscopic level, this swapping process generates the spatial gradient
of the temperature $T$ and the velocity $v_{x}$ along the $z$ axis in the
steady state. Moreover, we can calibrate the heat flux in the $z$-direction,
$q_{z}$, as

$\displaystyle{q}_{z}=\frac{\Delta E}{2tL_{x}L_{y}},$ (4)

where $\Delta E$ is the sum of the kinetic energy transferred by these
swapping processes during the simulation time $t$ and is calculated by the
following equation:

$\displaystyle\Delta
E=\sum_{t}\frac{m}{2}(v_{\text{h}}^{2}-v_{\text{c}}^{2}),$ (5)

where $v_{h}=\sqrt{v_{y}^{2}+v_{z}^{2}}$ is the kinetic energy of the particle
selected at the hot slab and $v_{c}=\sqrt{v_{y}^{2}+v_{z}^{2}}$ is that at the
cool slab in the swapping process. By using the heat flux $q_{z}$ calibrated
as Eq. (4) and the temperature gradient $\langle\frac{\partial T}{\partial
z}\rangle$ calibrated from the local temperature in each slab, the thermal
conductivity transverse to the shear flow is obtained from Eq. (1) as

$\displaystyle\kappa_{zz}=-\frac{q_{z}}{\langle\frac{\partial T}{\partial
z}\rangle}.$ (6)

In this work, MD simulations are performed using the LAMMPS package [20]. The
source code of the modified RNEMD method is available from the authors upon
request.

##  III SIMULATIONS

Simulations are performed for the simple LJ fluid and the KG model polymer
melt. In Sec. IIIA, we test the validity of our method by comparing our RNEMD
results with those obtained in a previous study[16]. In Sec. IIIB, we perform
RNEMD simulations for various lengths of polymer chains.

###  III.1 Lennard-Jones fluid

We consider the truncated LJ potential (the so-called Weeks-Chandler-Anderson
(WCA) potential [21]),

$\displaystyle
U_{\mathrm{LJ}}=\begin{cases}4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]+\epsilon&(r<2^{1/6}\sigma)\\\
0&(r\geq 2^{1/6}\sigma)\end{cases}\;.$ (7)

Hereafter, we measure the space and time in units of $\sigma$ and
$\sqrt{m\sigma^{2}/\epsilon}$, respectively. The temperature $T$ is measured
in units of $\epsilon/k_{\text{B}}$. The number density $\rho_{0}=0.85$ and
the sizes of the simulation cells $L_{x}=L_{y}=10.06$ and $L_{z}=40.24$ are
fixed. The rectangular simulation cell is divided into twenty slabs along the
$z$ axis. The equations of motion are integrated using the velocity Verlet
algorithm with a time step size of $\Delta t=0.001$.

The system is first equilibrated in the quiescent state at $T=1.25$ for
$10^{6}$ time steps with the Langevin thermostat. Next, swapping of the
kinetic energy is carried out with the time interval $\Delta
t_{\text{k}}=100\Delta t$ for the period $t=10^{6}\Delta t$ to achieve
stationary heat flux. Finally, the simulation involving heat and momentum
fluxes is performed for $10^{7}$ time steps in a variation of the time
interval of the momentum swapping as $\Delta t_{m}=$100$\Delta t$ –
1000$\Delta t$. The time interval of the kinetic energy swapping is fixed as
$\Delta t_{\text{k}}=100$$\Delta t$.

Figure 2:  The flow velocity profiles (a) and the temperature profiles (b) for
the LJ fluid are calculated with different intervals of momentum swapping
$\Delta t_{\text{m}}$, as shown in the legend to the figure. The interval of
the kinetic energy swapping is fixed at $\Delta t_{\text{k}}$=100$\Delta t$.
Figure 3:  Energy fluxes arising in momentum swapping (blue squares), energy
fluxes arising in kinetic energy swapping (orange circles), and total energy
fluxes (green triangles) are plotted as functions of the momentum swapping
interval. Each energy flux is measured in the lower half of the simulation
cell ($0<z<L_{z}/2$). Figure 4:  The thermal conductivity of the LJ fluid is
plotted as a function of the shear rate $\dot{\gamma}$. The blue squares show
the results obtained in this study, while the orange dashed line shows the
results obtained in Ref. [16] (i.e., Eq. (3)).

Figure 2 shows the spatial profiles of the flow velocity $U_{x}$ (in (a)) and
the temperature $T$ (in (b)) at different intervals of momentum swapping. The
flow velocity profiles are almost linear irrespective of the swapping interval
$\Delta t_{m}$. In contrast, the temperature profile becomes rather parabolic
as the frequency of momentum swapping $(\Delta t_{m})^{-1}$ increases for
$\Delta t_{m}\gtrsim 200$. This phenomenon was also reported in previous
studies[18, 22].

It is also observed that the temperature gradient decreases as the frequency
of momentum swapping increases. This behavior is explained in terms of the
contribution of momentum swapping to the energy flux. As mentioned earlier,
swapping atomic velocities $v_{x}$ between hot and cold slabs will likely
generate a negative energy flux along the $z$ axis in the lower half of the
simulation cell ($0<z<L_{z}/2$). In Fig. 3, we show the energy fluxes due to
the swapping of kinetic energy and those due to the swapping of momentum
separately. It is clear that the energy flux due to kinetic energy swapping is
almost constant for $\Delta t_{m}\gtrsim 200$$\Delta t$. In contrast, the
energy flux due to momentum swapping, which has a negative value, becomes
increasingly prominent as the swapping frequency $(\Delta t_{m})^{-1}$
increases when $\Delta t_{m}\lesssim 500$$\Delta t$. Thus, the net energy flux
decreases as the frequency of momentum swapping increases, causing the
temperature gradient to decrease.

Figure 4 shows the dependency of the shear rate $\dot{\gamma}$ on the thermal
conductivity in the direction transverse to shear flows. In this figure, we
plot the deviation of $\kappa_{zz}$ from that in the quiescent state,
$\kappa_{zz}-\kappa_{0}$, where $\kappa_{0}$ is the thermal conductivity
obtained in the quiescent state as a function of the shear rate
$\dot{\gamma}$. Our results are in good agreement with the values obtained in
a previous study [16] and consistent with the theoretically expected quadratic
$\dot{\gamma}$ dependence [15] within the wide range of shear rates, even
though at small swapping intervals (e.g., $\Delta t_{m}$=100 and 120), the
temperature profiles become rather parabolic. This result demonstrates that
our RNEMD method is valid within the wide range of $\dot{\gamma}$ shown in
Fig. 4. Here, we also emphasize that $\kappa_{zz}$ remains almost constant at
low shear rates (e.g., $\dot{\gamma}<0.025$) in the LJ fluid.

###  III.2 Polymer melt

We consider the KG model polymer melt. The simulation system includes
$N_{\text{c}}$ polymer chains with $N_{\text{b}}$ beads. All bead practices
interact via the repulsive part of the LJ potential (7), and the consecutive
beads on each chain are connected by an anharmonic spring potential

$\displaystyle
U_{\mathrm{FENE}}=-\frac{1}{2}k_{c}R^{2}_{0}\ln\left[1-\left(\frac{r}{R_{0}}\right)^{2}\right],$
(8)

where $k_{\text{c}}=30\epsilon/\sigma^{2}$ and $R_{0}=1.5\sigma$ are fixed.
The number density $\rho_{0}=0.85$ and the initial temperature $T_{0}=1$ are
fixed. We carried out simulations of model polymer melts with different chain
lengths $N_{\text{b}}$. In Table I, the parameter sets used in the simulations
are summarized.

Table 1: Parameter values used in the simulations for KG model polymer melts are shown. The first and second columns show the number of beads per molecule and the number of molecules, respectively. The third to fifth columns show the lengths of the simulation cell in the x-direction, y-direction, and z-direction, respectively.  $N_{\text{b}}$ | $N_{c}$ | $L_{x}$ | $L_{y}$ | $L_{z}$
---|---|---|---|---
30 | 1600 | 24.16 | 24.16 | 96.64
50 | 1100 | 25.29 | 25.29 | 101.16
100 | 600 | 26.03 | 26.03 | 104.12
200 | 360 | 27.66 | 27.66 | 110.64

Each simulation is performed for $10^{8}$ time steps with $\Delta t=0.001$.
The time interval of momentum swapping $\Delta t_{m}$ is between 100 and 600
time steps, while the time interval of kinetic energy swapping $\Delta
t_{\text{k}}=200\Delta t$ is fixed. Within this range of $\Delta t_{m}$, the
shear rates generated in the RNEMD simulations are at most less than 0.005. In
the present study, for polymer melts, we only focus on the linear regime where
the higher-order term of the shear rate is neglected.

Figure 5 shows the spatial profiles of flow velocity $U_{x}$ and temperature
$T$ for polymer melts with different chain lengths. Both the velocity and
temperature profiles are nearly linear within the parameter range.

Figure 5:  The flow velocity profiles (a) and the temperature profiles (b) calculated for the KG model polymer melt with different momentum swapping intervals $\Delta t_{\text{m}}$=100$\Delta t$ — 600$\Delta t$. The KE swapping interval is fixed at $\Delta t_{\text{k}}$=200$\Delta t$. Figure 6:  The thermal conductivity plotted as the function of shear rate $\dot{\gamma}$ for the KG model polymer melts. The red pentagons show the results for $N_{\text{b}}=30$, the blue squares for $N_{\text{b}}=50$, the green points for $N_{\text{b}}=100$, and the orange triangles for $N_{\text{b}}=200$. Table 2: Thermal conductivities and normal stresses in the quiescent state for the KG model polymer melts with different lengths of polymer chains. The first column is the number of beads per molecule. The second and the third columns are thermal conductivity and normal stress, respectively. $N_{\text{b}}$ |  $\kappa_{0}$ |  $\tau_{0}$
---|---|---
30 |  4.97 |  -5.00
50 |  4.94 |  -4.91
100 |  4.93 |  -4.89
200 |  4.93 |  -4.89

Figure 6 shows the thermal conductivity $\kappa_{zz}$ as a function of the
shear rate $\dot{\gamma}$. Here, $\kappa_{0}$ is the thermal conductivity in
the quiescent state, and the values of $\kappa_{0}$ for different chain
lengths are summarized in Table 2. For $N_{\text{b}}=$30 and
$N_{\text{b}}=$50, the thermal conductivity remains almost unchanged even with
an increase in $\dot{\gamma}$. However, for $N_{\text{b}}=$100 and
$N_{\text{b}}=$200, the thermal conductivity clearly decreases as the shear
rate increases.

Figure 7 shows the relationship between the normal stress difference and bond
orientation (in the upper figure) and the relationship between the thermal
conductivity and bond orientation (in the lower figure). In each figure, the
dashed lines represent the linear fits of the simulation results for each
polymer chain. Here, the bond-orientation tensor $Q_{ij}$ ($i,j=\\{x,y,z\\}$)
is defined as

$\displaystyle
Q_{ij}=\frac{1}{N_{\text{c}}}\sum_{\text{chains}}\frac{1}{N_{\text{b}}-1}\sum_{k=1}^{N_{\text{b}}-1}\frac{b_{i}^{k}}{b_{\text{min}}}\frac{b_{j}^{k}}{b_{\text{min}}},$
(9)

where ${b_{i}^{k}}$ for $1\leq k\leq N_{\text{p}}-1$ is the component of the
bond vector between consecutive beads in the same chain and $b_{\text{min}}$
is the distance at which the total potential
$U_{\text{LJ}}(r)+U_{\text{FENE}}(r)$ is minimized (i.e.,
$b_{\text{min}}\simeq 0.97$). On the horizontal axis, $Q_{0}$ is the value in
the quasi-simulation system (i.e., $Q_{0}\simeq 0.33$). The stress tensor
$\tau_{ij}$ ($i,j=\\{x,y,z\\}$) is defined as

$\displaystyle\tau_{ij}=-\frac{1}{V}\sum_{k}\left[m(v_{i}^{k}-U_{i}^{k})(v_{j}^{k}-U_{j}^{k})+\frac{1}{2}\sum_{l\neq
k}F_{i}^{kl}r_{j}^{kl}\right],$ (10)

where $v_{i}^{k}$ and $U_{i}^{k}$ denote the atomic velocity of the $k$th bead
and the macroscopic flow velocity at the position of the $k$th bead,
respectively. The force acting on the $k$th bead by the $l$th bead is written
as $F_{i}^{kl}$. Here, $\tau_{0}$ is the normal stress in the quiescent state,
and the values of $\tau_{0}$ for polymer melts with different chain lengths
are summarized in Table 2.

The normal stress difference is linearly proportional to the bond orientation
for all the polymer chains with different numbers of beads $N_{\text{b}}$.
Thus, we obtain the well-known stress-optical rule given by

$\displaystyle\frac{\tau_{0}-\tau_{zz}}{|\tau_{0}|}=C_{\text{s}}(Q_{0}-Q_{zz}),$
(11)

where $C_{\text{s}}$ is the stress-optical coefficient. The values of
$C_{\text{s}}$ are calculated as 0.32, 0.39, 0.39, and 0.37 for
$N_{\text{b}}=30$ $N_{\text{b}}=50$, $N_{\text{b}}=100$ and
$N_{\text{b}}=200$, respectively. Our simulation results show that the linear
stress-optical rule robustly holds irrespective of the length of the polymer
chain even if a small thermal gradient exists.

Figure 7:  The deviations in the normal stress (a) and thermal conductivity
(b) from those in the quiescent state are plotted as functions of the
deviation of the bond orientation from that in the quiescent state. The red
squares show the results for the KG model polymer melt with $N_{\text{b}}=30$,
the blue squares for $N_{\text{b}}=50$, the green downward triangles for
$N_{\text{b}}=100$, and the orange upward triangles for $N_{\text{b}}=200$.

In contrast, for the relationship between the thermal conductivity and the
bond orientation, the linear dependency of the bond orientation on the thermal
conductivity is not clearly observed for short chains $N_{\text{b}}=30$ and
$N_{\text{b}}=50$, whereas it appears that the linear relationship holds for
long chains $N_{\text{b}}$=100 and $N_{\text{b}}=$200, as expressed by

$\displaystyle\frac{\kappa_{0}-\kappa_{zz}}{\kappa_{0}}=C_{\text{b}}(Q_{0}-Q_{zz}),$
(12)

where $C_{\text{b}}$ is the thermal-optical coefficient. Furthermore, once the
linear relationship holds at $N_{\text{b}}\gtrsim 100$, $C_{\text{b}}$ is not
affected by the beard number $N_{\textrm{b}}$. The values of $C_{\text{b}}$
are calculated as 0.67 and 0.71 for $N_{\text{b}}=100$ and $N_{\text{b}}=200$,
respectively.

Notably, the present KG model polymer melt represents unentangled polymers at
$N_{\text{b}}$=30 and 50 but entangled polymers at $N_{\text{b}}$=100 and
$N_{\text{b}}=$200 (see Appendix). This indicates that the linear relationship
between the thermal conductivity and the bond-orientation tensor holds only
for entangled polymers and does not hold for unentangled polymers.

This behavior is intuitively explained as follows: Microscopically, the
thermal conductivity of a polymer melt is determined mainly via two
contributions: the interaction of inconsecutive bead particles (as in the LJ
fluid) and the transmission of molecular vibrations in consecutive bead
particles along polymer chains. The results of the LJ fluid in Sec. IIIA
indicate that the former contribution does not affect the thermal conductivity
for the small strain rates considered in this section. Thus, it is thought
that the deviation of the thermal conductivity from that in the quiescent
state, $\kappa_{0}-\kappa_{zz}$, stems from the latter contribution. Once
entanglement is generated between polymer chains, a network of polymer chains
is constructed, and the transmission of molecular vibrations is more likely to
occur along the network structure. Since the network structure is
characterized by the bond orientation tensor, as in the STR derived from
network theory [3], the deviation of the thermal conductivity
$\kappa_{0}-\kappa_{zz}$ linearly depends on that of the bond orientation
tensor $Q_{0}-Q_{zz}$ when entanglement is generated between polymer chains.

Finally, Figure 8 shows the relation between the thermal conductivity
$\kappa_{0}-\kappa_{zz}$ and the stress difference $\tau_{0}-\tau_{zz}$
expressed by

$\displaystyle\frac{\kappa_{0}-\kappa_{zz}}{\kappa_{0}}=C_{\text{t}}G_{\text{N}}\frac{\tau_{0}-\tau_{zz}}{G_{\text{N}}},$
(13)

where $G_{\text{N}}$ is the plateau modulus and is estimated to be $8.0\times
10^{-3}$ for our simulation model (see Appendix). The above equation is
slightly different from the STR described by Eq. (2). In Eq. (13), the trace
of the thermal conductivity tensor $\mathrm{tr}(\bm{\kappa})$ in Eq. (2) is
replaced with $\kappa_{0}$.

Figure 8: The deviation of the thermal conductivity from that in the quiescent
state is plotted as a function of the deviation of the normal stress from that
in the quiescent state for entangled model polymer melts. The green downward
triangles show the results for $N_{\text{b}}=100$, and the orange upper
triangles show the results for $N_{\text{b}}=200$.

For entangled polymer chains $N_{\text{b}}=$100 and $N_{\text{b}}=$200, the
linear relationship (13) holds, and their coefficients
$C_{\text{t}}G_{\text{N}}$ are almost the same. This behavior is consistent
with the experimental results obtained in Refs. [23, 7, 6, 8, 24], where the
STR is measured under deformation in quasi-quiescent states after the
cessation of external forces. However, in our results, the value of the
coefficient is an order of magnitude smaller than the experimental result
$C_{\text{t}}G_{\text{N}}\simeq 0.04$ obtained in the references. In our
results, the values of the coefficient are calculated as
$C_{\text{t}}G_{\text{N}}=3.1\times 10^{-3}$ and $3.3\times 10^{-3}$ for
$N_{\text{b}}=100$ and $N_{\text{b}}=200$, respectively. The reason for this
discrepancy has not been explained.

##  IV CONCLUDING REMARKS

We investigated the STR of polymer melts in the transverse direction under
shear flows by using molecular dynamics simulations. We extended the original
RNEMD method to simultaneously generate spatial gradients of flow velocity and
temperature in a single direction. This method enables the accurate
measurement of thermal conductivity in the transverse direction to shear flow.
The modified RNEMD method was first verified through comparisons with
theoretical and numerical results, given in the references, on the nonlinear
dependence of the thermal conductivity on the shear rate for LJ fluids.
Subsequently, this approach was applied to KG model polymer melts.

Our main finding is that the STR holds for entangled polymer melts even under
shear flows, but it does not hold for unentangled polymer melts. This behavior
completely differs from the stress-optical rule (which uniformly holds for
polymer melts irrespective of chain length). Our observation that the
occurrence of an STR is linked to the occurrence of polymer chain entanglement
aligns with the foundation of an STR since an STR was originally derived based
on the network theory of polymer melts that considers energy transmission
along the network structure of entangled polymer chains.

Furthermore, our numerical results demonstrated that once entanglement occurs
in polymer chains, the stress-thermal coefficient is not affected by the chain
length. This observation is consistent with the experimental results of
anisotropic thermal conductivity measured in a quasi-quiescent state following
shear or elongational deformations. However, our numerical results indicate
that the stress-thermal coefficient $C_{t}G_{N}$ in Eq. (13) is one order of
magnitude smaller in shear flows than in the quasi-quiescent state. The
mechanism of the decrease in the stress-thermal coefficient in shear flows has
yet to be explained. However, our observation seems reasonable because the
energy transmission along the network of polymer chains should be hindered by
successive breakages of the network structure during shear flow.

Apart from the physical aspects of the STR, this study fundamentally
contributes to the computer simulation of the thermo-fluid dynamics of polymer
melts. Instead of the conventional Fourier’s law with an isotropic thermal
conductivity, by utilizing either the STR or the thermal-optical relation as
the constitutive relation for the energy flux, one can compute the coupled
dynamics of flow behavior and temperature distribution more accurately. This
aspect is highly important for computational thermo-fluid dynamics of polymer
melts with board engineering applications, for example, polymer processing in
chemical engineering and lubrication systems in mechanical engineering. In a
future study, the distinct behaviors produced by the STR in thermo-fluid
dynamics will be addressed.

Finally, in the present study, we considered only the STR in the transverse
direction to the shear flow. However, the original STR is given in tensor
form. Thus, to confirm the STR, we need to investigate the longitudinal
direction of the shear flow. This will also be addressed in our future study.

*

##  Appendix A Relaxation functions of unentangled and entangled polymer melts

In the appendix, we investigate the differences in relaxation behavior between
the model polymer melts composed of short chains (i.e., $N_{\text{b}}=30$ and
50) and those composed of long chains (i.e., $N_{\text{b}}=100$ and 200) in
the equilibrium state. Furthermore, we present the method for estimating
$G_{\text{N}}$ in Eq. (13). The relaxation dynamics of short-chain models with
$N_{\text{b}}<N_{\text{e}}$, where $N_{\text{e}}$ represents the average bead
number between consecutive entanglements, can be described by Rouse dynamics,
while those of long-chain models with $N_{\text{b}}>N_{\text{e}}$ can be
described by reptation dynamics [25]. The value of $N_{\text{e}}$ is typically
in the range of 60 to 100 [19, 26, 27, 28].

Figure 9 shows the normalized time-correlation function of the end-to-end
vector of the polymer chain $C(t)$ calculated as

$\displaystyle
C(t)=\frac{\langle\mathbf{P}(t+t_{0})\cdot\mathbf{P}(t_{0})\rangle}{\langle\
|\mathbf{P}_{(}t_{0})|^{2}\rangle},$ (14)

where $\mathbf{P}(t)$ represents the end-to-end vector of each polymer chain.
Figure 10 shows the stress relaxation function $G(t)$ calculated as

$\displaystyle
G(t)=\frac{V}{3k_{\text{B}}T}(\langle\tau_{xy}(t+t_{0})\tau_{xy}(t_{0})\rangle+\langle\tau_{xz}(t+t_{0})\tau_{xz}(t_{0})\rangle+\langle\tau_{yz}(t+t_{0})\tau_{yz}(t_{0})\rangle).$
(15)

According to $C(t)$, both the Rouse dynamics and the reptation dynamics
predict the same simple functional form in terms of the relaxation time
$\tau_{d}$ as

$\displaystyle
C(t)=\sum_{\text{odd}\>p}\frac{8}{\pi^{2}p^{2}}\exp(-p^{2}t/\tau_{d}),$ (16)

where the summation is over odd $p$ and $1\leq p\leq N_{\text{b}}-1$. Since
the first term $p=1$ in the summation is dominant for any $t$, this function
decays nearly exponentially, and $\tau_{\text{d}}$ may be estimated from
$C(\tau_{d})\simeq\exp(-1)$. We obtained the relaxation times $\tau_{d}$ for
$N_{\text{b}}=30,50,100$ and $200$ as $\tau_{d}$=740, 2380, 13240, and 73980.

Theoretically, in Rouse dynamics, the stress relaxation function
$G_{\text{Rouse}}(t)$ is expressed as of the relaxation time $\tau_{\text{R}}$
as

$\displaystyle G_{\text{Rouse}}(t)=\frac{\rho
T}{N_{\text{b}}}\sum_{p=1}^{N_{\text{b}}-1}\exp(-p^{2}t/\tau_{R}N_{\text{b}}^{2})$
(17)

and in reptation dynamics, the stress relaxation function $G_{\text{rep}}(t)$
is given by [25]

$\displaystyle G_{\text{rep}}(t)=\frac{4\rho T}{5N_{\text{e}}}C(t).$ (18)

Here, the models with $N_{\text{b}}=30$ and $N_{\text{b}}=50$ are assumed to
be represented by Rouse dynamics. The relaxation time $\tau_{d}$ for
$N_{\text{b}}=30$ should be equal to the Rouse relaxation time as follows

$\displaystyle\tau_{\text{d}}=\tau_{R}N_{\text{b}}^{2}.$ (19)

We obtained $\tau_{R}=0.82$ for $N_{\text{b}}=30$ and $N_{\text{b}}=50$,
assuming that $\tau_{R}$ does not depend on $N_{\text{b}}$.

On the other hand, when we suppose that the relaxation dynamics of models with
$N_{\text{b}}=100$ and $N_{\text{b}}=200$ are described by the reptation
dynamics, $N_{\text{e}}$ is estimated as $N_{\text{e}}=85$ for both
$N_{\text{b}}=100$ and $N_{\text{b}}=200$ by fitting Eq. (18) to our numerical
results in Fig. 10 at $t=\tau_{d}$. In Figure (10), we also compare our
simulation results and the theoretical formulas with the fitting values
obtained above. Our numerical results $G(t)$ for $N_{\text{b}}=30$ and
$N_{\text{b}}=50$ exhibit Rouse dynamics, while those for $N_{\text{b}}=100$
and $N_{\text{b}}=200$ closely match the reptation dynamics. The plateau
modulus $G_{\text{N}}$ is determined by the following equation:

$\displaystyle G_{\text{N}}=\frac{4\rho T}{5N_{\text{e}}}.$ (20)

Thus, the value of $G_{\text{N}}$ is estimated as $8.0\times 10^{-3}$ for both
$N_{\text{b}}=100$ and $N_{\text{b}}=200$.

Figure 9:  The normalized time-correlation function of the end-to-end vectors
$C(t)$ for the KG model polymer melts with $N_{\text{b}}=30$ (red line),
$N_{\text{b}}=50$ (blue line), $N_{\text{b}}=100$ (green line) and
$N_{\text{b}}=200$ (orange line).  Figure 10: The stress autocorrelation
function $G(t)$ for the KG model polymer melts with $N_{\text{b}}=30$ (red
solid line), $N_{\text{b}}=50$ (blue solid line), $N_{\text{b}}=100$ (green
solid line) and $N_{\text{b}}=200$ (orange solid line). The Rouse relaxation
function $G_{\text{Rouse}}(t)$ in Eq. (17) for $N_{\text{b}}=30$ (red dotted
line) and $N_{\text{b}}=50$ (blue dotted line) and the reptation relaxation
function $G_{\text{rep}}(t)$ in Eq. (18) for $N_{\text{b}}=100$ (green dashed
line) and $N_{\text{b}}=200$ (orange dashed line) are also plotted.

###### Acknowledgements.

The computations in this work were performed using the supercomputers of
ACCMS, Kyoto University and the Center for Cooperative Work on Data Science
and Computational Science, University of Hyogo. This work was supported by
JSPS KAKENHI Grant Number 23H03350.

## References

  * Yasuda and Yamamoto [2014] S. Yasuda and R. Yamamoto, Synchronized molecular-dynamics simulation via macroscopic heat and momentum transfer: An application to polymer lubrication, Physical Review X 4, 041011 (2014).
  * Yasuda [2019] S. Yasuda, Synchronized molecular-dynamics simulation of the thermal lubrication of an entangled polymeric liquid, Polymers 11, 131 (2019).
  * Van Den Brule [1989] B. H. A. A. Van Den Brule, A network theory for the thermal conductivity of an amorphous polymeric material, Rheologica Acta 28, 257 (1989).
  * loong Choy  _et al._ [1978] C. loong Choy, W. H. Luk, and F. C. Chen, Thermal conductivity of highly oriented polyethylene, Polymer 19, 155 (1978).
  * loong Choy  _et al._ [1999] C. loong Choy, Y. W. Wong, G. W. Yang, and T. Kanamoto, Elastic modulus and thermal conductivity of ultradrawn polyethylene, Journal of Polymer Science Part B 37, 3359 (1999).
  * Simavilla  _et al._ [2012] D. N. Simavilla, J. D. Schieber, and D. C. Venerus, Anisotropic Thermal Transport in a Crosslinked Polyisoprene Rubber Subjected to Uniaxial Elongation, Journal of Polymer Science Part B: Polymer Physics , 1 (2012).
  * Schieber  _et al._ [2012] J. D. Schieber, D. C. Venerus, and S. Gupta, Molecular origins of anisotropy in the thermal conductivity of deformed polymer melts: stress versus orientation contributions, Soft Matter 8, 11781 (2012).
  * Simavilla and Venerus [2014] D. N. Simavilla and D. C. Venerus, Investigation of anisotropic thermal conductivity in polymers using infrared thermography, Journal of Heat Transfer 136, 10.1115/1.4028324 (2014).
  * Broerman  _et al._ [1999] A. W. Broerman, D. C. Venerus, and J. D. Schieber, Evidence for the stress-thermal rule in an elastomer subjected to simple elongation, Journal of Chemical Physics 111, 6965 (1999).
  * Nieto Simavilla  _et al._ [2020] D. Nieto Simavilla, A. P. Sgouros, G. G. Vogiatzis, C. Tzoumanekas, V. Georgilas, W. M. Verbeeten, and D. N. Theodorou, Molecular dynamics test of the stress-thermal rule in polyethylene and polystyrene entangled melts, Macromolecules 53, 789 (2020).
  * Iddir  _et al._ [2000] H. Iddir, D. C. Venerus, and J. D. Schieber, Measuring anisotropic thermal conduction in polyisobutylene following step shear strains, AIChE Journal 46, 610 (2000).
  * Schieber  _et al._ [2004] J. D. Schieber, D. C. Venerus, K. Bush, V. Balasubramanian, and S. Smoukov, Measurement of anisotropic transport in flowing polymers by using a holographic technique, Proceedings of the National Academy of Sciences of the United States of America 101, 13142 (2004).
  * Venerus  _et al._ [2004] D. C. Venerus, J. D. Schieber, V. Balasubramanian, K. Bush, and S. Smoukov, Anisotropic thermal conduction in a polymer liquid subjected to shear flow, Phys. Rev. Lett. 93, 098301 (2004).
  * Dai and Tanner [2006] S. C. Dai and R. I. Tanner, Anisotropic thermal conductivity in sheared polypropylene, Rheologica acta 45, 228 (2006).
  * Daivis and Evans [1993] P. J. Daivis and D. J. Evans, Thermal conductivity of a shearing fluid, Phys. Rev. E 48, 1058 (1993).
  * Smith  _et al._ [2019] E. R. Smith, P. J. Daivis, and B. D. Todd, Measuring heat flux beyond Fourier law, The Journal of Chemical Physics 150 (2019).
  * Müller-Plathe [1999] F. Müller-Plathe, Reversing the perturbation in nonequilibrium molecular dynamics: An easy way to calculate the shear viscosity of fluids, Physical Review E 59, 4894 (1999).
  * Müller-Plathe [1997] F. Müller-Plathe, A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity, The Journal of chemical physics 106, 6082 (1997).
  * Kremer and Grest [1990] K. Kremer and G. S. Grest, Dynamics of entangled linear polymer melts: A molecular-dynamics simulation, The Journal of Chemical Physics 92, 5057 (1990).
  * Plimpton [1995] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117, 1395 (1995).
  * Weeks  _et al._ [2003] J. D. Weeks, D. Chandler, and H. C. Andersen, Role of Repulsive Forces in Determining the Equilibrium Structure of Simple Liquids, The Journal of Chemical Physics 54, 5237 (2003).
  * Tenney and Maginn [2010] C. M. Tenney and E. J. Maginn, Limitations and recommendations for the calculation of shear viscosity using reverse nonequilibrium molecular dynamics, J. Chem. Phys. 132, 1395 (2010).
  * Gupta  _et al._ [2013] S. Gupta, J. D. Schieber, and D. C. Venerus, Anisotropic thermal conduction in polymer melts in uniaxial elongation flows, Journal of Rheology 57, 427 (2013).
  * Venerus  _et al._ [2001] D. C. Venerus, J. D. Schieber, I. H., G. J., and B. A., Anisotropic thermal diffusivity measurements in deforming polymers and the stress-thermal rule, International Journal of Thermophysics 22, 1215 (2001).
  * Doi and Edwards [1989] M. Doi and S. F. Edwards, The theory of polymer dynamics, Journal of Polymer Science Part C: Polymer Letters 27, 239 (1989).
  * Lee and Kremer [2009] W. B. Lee and K. Kremer, Entangled polymer melts: Relation between plateau modulus and stress autocorrelation function, Macromolecules 42, 6270 (2009).
  * Pütz  _et al._ [2000] M. Pütz, K. Kremer, and G. S. Grest, What is the entanglement length in a polymer melt?, Europhysics Letters 49, 735 (2000).
  * Yamamoto and Onuki [2004] R. Yamamoto and A. Onuki, Entanglements in quiescent and sheared polymer melts, Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 70, 12 (2004).