Solidification of the Lennard-Jones fluid near the wall in thermohydrodynamic lubrication

TL;DR

This study uses molecular dynamics simulations to reveal that viscous heating can induce solidification of Lennard-Jones fluids near walls in thermohydrodynamic lubrication, highlighting a counterintuitive phase transition.

Contribution

It uncovers the phenomenon of viscous heating-induced solidification of LJ fluids near walls and analyzes its microscopic and macroscopic mechanisms.

Findings

Solidification occurs only with sufficient viscous heating.

LJ molecules densely confine near the wall and solidify upon interaction with crystallized wall molecules.

Band formation appears in highly confined regimes.

Abstract

We investigate the thermohydrodynamic lubrication of the Lennard-Jones (LJ) fluid in the parallel-plate channel composed of the LJ particles by using molecular dynamics (MD) simulation. We discover a counterintuitive solidification of the LJ fluid near the wall, i.e., \textit{viscous heating-induced solidification}, where solidification occurs only when the viscous heating of the LJ fluid is sufficiently large. The solidification mechanism is investigated from both macroscopic and microscopic points of view. It is found that the LJ molecules are densely confined in the vicinity of the wall via the thermohydrodynamic transport of the bulk fluid and that when the local density in the vicinity of the wall is close to the solidification line in the phase diagram, the LJ molecules are solidified due to the interaction with the crystallized wall molecules. Band formation is also observed in the highly confined regime when the channel width is sufficiently large.

Figures

Equations

$$U(r)=\left\{\begin{array}[]{cc}4\varepsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right],&(0<r<r_{c}),\\ 0,&(r_{c}\leq r).\end{array}\right.$$

$$m{\ddot{\mathbf{r}}}_i(t)=-\sum _j\frac{\partial U(|{\mathbf{r}}_{ij}|)}{\partial {\mathbf{r}}_{ij}}$$

$$m{\ddot{\mathbf{r}}}_i(t)=-\sum _j\frac{\partial U(|{\mathbf{r}}_{ij}|)}{\partial {\mathbf{r}}_{ij}}-k_s({\mathbf{r}}_i-{\mathbf{r}}_i^w)-\gamma \dot{{\mathbf{r}}_i}+\mathbf{R}(t),$$

$$<R_{\alpha }(t)R_{\beta }(t-s)>=2m{k}_B{T}^w\gamma {\delta }_{\alpha \beta }\delta (s).$$

$$\rho (l)=\frac{1}{|V_{\mathrm{bin}}^l|}{\int }_{\mathbf{r}\in V_{\mathrm{bin}}^l}\sum _i\delta (\mathbf{r}-{\mathbf{r}}_i)d\mathbf{r},$$

$$\rho (l)u_x(l)=\frac{1}{|V_{\mathrm{bin}}^l|}{\int }_{\mathbf{r}\in V_{\mathrm{bin}}^l}\sum _i{\dot{r}}_{xi}\delta (\mathbf{r}-{\mathbf{r}}_i)d\mathbf{r},$$

$$\rho (l)T(l)=\frac{1}{3|V_{\mathrm{bin}}^l|}{\int }_{\mathbf{r}\in V_{\mathrm{bin}}^l}\sum _i\dot{(}{\mathbf{r}}_i-u_x(l){\delta }_{x\alpha }{)}^2\delta (\mathbf{r}-{\mathbf{r}}_i)d\mathbf{r},$$

$$p_{\alpha \beta }(l)=\frac{1}{|V_{\mathrm{bin}}^l|}{\int }_{\mathbf{r}\in V_{\mathrm{bin}}^l}\sum _i\left(({\dot{r}}_{\alpha i}-u_x(l){\delta }_{x\alpha }){\dot{r}}_{\beta i}+r_{\alpha i}{f}_{\beta i}\right)\delta (\mathbf{r}-{\mathbf{r}}_i)d\mathbf{r},$$

$$g^l(r)={⟨\frac{n_i(r)}{4\pi r^{2}dr{\rho }_0}⟩}_l,$$

$$D_y={\int }_0^{\infty }{⟨\overline{v_y^i(t+{\tau }_0)v_y^i({\tau }_0)}⟩}_{l}dt,$$

$$\frac{\mu V_w^2/H^2}{\lambda \Delta T/H^2}\propto V_w^2,$$

$$-\frac{d^{2}T}{d{\hat{y}}^2}=\frac{\mu }{\lambda }{V}_w^2,$$

$$T(\hat{y})=4\Delta T\hat{y}(1-\hat{y})+T_w,$$

$${\int }_0^1\rho (\hat{y})d\hat{y}={\rho }_0,$$

$${\int }_0^1\mathcal{F}(p,T(\hat{y}))d\hat{y}={\rho }_0.$$

$$\begin{split}&\int_{0}^{1}\left(\frac{dp}{dV_{w}}\right)\left(\frac{\partial{\cal F}}{\partial p}\right)_{T}d\hat{y}+\int_{0}^{1}\left(\frac{\partial T}{\partial V_{w}}\right)\left(\frac{\partial{\cal F}}{\partial T}\right)_{T}d\hat{y}=0,\\ &\frac{dp}{dV_{w}}=-\int_{0}^{1}\left(\frac{\partial T}{\partial V_{w}}\right)\left(\frac{\partial{\cal F}}{\partial T}\right)_{T}d\hat{y}\Big{/}\int_{0}^{1}\left(\frac{\partial{\cal F}}{\partial p}\right)_{T}d\hat{y}\quad>0.\end{split}$$

$$\rho (l)=\frac{1+{\int }_0^{r_c}4\pi r^2{\rho }_0{g}^l(r)dr}{\frac{4}{3}\pi r_c^3},$$

Full text

I Introduction

Lubrication in high-speed mechanical systems may involve complicated phenomena such as thermorheological coupling in the bulk fluid and chemical reactions and phase transitions at the interface. Computer simulation of such a complicated lubrication system is challenging and important from both engineering and scientific points of view.

For hydrodynamic lubrication problems, the computational fluid dynamics (CFD) approach is usually utilized by employing any constitutive relations, which stem from the microscopic dynamics of the molecules. Since constitutive relations are usually unknown for complex fluids involved in, e.g., thermorheological coupling, chemical reactions, and phase transitions, it is difficult to apply CFD simulation to complicated lubrication problems ab initio. Molecular dynamics (MD) simulation is applicable even for complex fluids because the complicated transport phenomena are autonomously reproduced once their molecular models are appropriately specified.book:89AT ; book:08EM ; book:05KBA

In this study, we analyze the high-speed lubrication of the Lennard-Jones (LJ) fluid in relation to thermohydrodynamic coupling and phase transitions. Instead of conventional fluid simulations, we perform MD simulation of thermal flows in a channel made of molecularly constituted thermal walls. Thus, thermohydrodynamic coupling, e.g., viscous heating, and phase transition are autonomously reproduced, as naturally occurs in real high-speed lubrication systems.

It is well known that simple liquids may be solidified in molecularly thin layers due to the confinement, i.e., confinement-induced phase transition art:95KK ; art:98KK ; art:06Aetal . Recently, we reported in a proceedings paper that a counterintuitive phase transition, i.e., viscous heating-induced solidification, occurs at the interface in the high-speed lubrication of the LJ fluid, where the channel width is much larger than the molecular size. art:18YO . In this paper, we comprehensively carry out MD simulations with changing channel width, wall speed, wall temperature, fluid density, and wall structure, and aim to unveil the mechanism of the solidification from both microscopic and macroscopic points of view.

MD simulations of thermohydrodynamic lubrication in nanochannels were previously reported in the literatures art:97KPY ; art:06HO ; art:10KBC ; art:17GB . These studies clarified the distinctive features of nanoscale flows, such as the slip on the boundaries. In this study, we consider a channel that is much larger than the molecular size, at which macroscopic transport is significant, and focus on the thermohydrodynamic coupling in the bulk fluid and the phase transition at the boundaries. The target of our study is more relevant to micromechanical engineering involving high-speed mechanical systems.

Incidentally, the multiscale hybrid method of MD and CFD is currently a very active research field art:03EE ; art:03KGHKRT ; art:05RE ; art:09KS ; art:13BLR ; art:16ZRE ; art:08YY ; art:09YY ; art:10YY ; art:11YY ; art:13MYTY . The thermohydrodynamic lubrication of polymeric fluids was also investigated by the synchronized MD method, in which a transitional behavior of the polymer conformation due to thermohydrodynamic coupling was clarified art:14YY ; art:16YY ; art:19Y . Although the multiscale approach is powerful and promising for complex fluids, to further develop the multiscale method, we require more first-principles results obtained through MD simulations. This study also aims to contribute to this objective.

In the following text, we describe the problem and simulation model considered in this study in Sec. II. The simulation results are given in Sec. III, where the counterintuitive solidification is revealed. In Sec. IV, the solidification mechanism is discussed from both macroscopic and microscopic points of view. Finally, we give the concluding remarks in Sec. V.

II Problem and simulation model

The LJ fluid between parallel plates, as shown in Fig. 1(a), is considered. The fluid domain extends $0<y<H$, and the wall domains extend $-W\leq y\leq 0$ (the lower wall) and $H\leq y\leq H+W$ (the upper wall).

Both the fluid and the walls are composed of LJ particles that interact with each other via the LJ potential:

[TABLE]

Here, $r_{c}$ is the cut-off parameter, and $\varepsilon$ and $\sigma$ are the units of energy and length of the LJ particles, respectively.

The wall particles are connected to the face-centered cubic (FCC) lattice structure $\{\bm{r}^{w}_{i}\}$ by the spring potential, and the temperature of the wall particles is kept at a constant value of $T^{w}$ by the Langevin thermostat algorithm.

Thus, the dynamics of the LJ particles are described, for the fluid particles (i.e., ${r_{y}}_{i}\in(0,H)$), by

[TABLE]

and, for the wall particles (i.e., ${r_{y}}_{i}\in[-W,0]\cup[H,H+W]$), by

[TABLE]

where $R_{\alpha}(t)$ ($\alpha=x,y,z$) is the white Gaussian noise that satisfies

[TABLE]

Here, $\bm{r}_{i}$ represents the position of the $i$th particle, $\bm{r}_{ij}$ is defined as $\bm{r}_{ij}=\bm{r}_{i}-\bm{r}_{j}$, $m$ is the mass of an LJ particle, $k_{s}$ is the spring constant, $\gamma$ is the damping coefficient, and $k_{B}$ is the Boltzmann constant. We note that the summation $\sum_{j}$ applies to both the fluid and wall particles. We also note that the wall and fluid particles are considered the same in size and mass.

The cut-off length $r_{c}=2.8$, the spring constant $k_{s}=10$ and the damping coefficient $\gamma=0.1$ are fixed. Hereafter, we express quantities in units of mass $m$, energy $\varepsilon$, length $\sigma$, and time $\tau=\sqrt{m\sigma^{2}/\varepsilon}$.

We note that the temperature of the fluid domain is not artificially controlled by any thermostat algorithm but varies autonomously according to the mass, momentum, and heat transfer between the parallel plates. Only the walls are kept at a constant temperature $T^{w}$ by the Langevin thermostat as Eq. (3).

The LJ fluid is initially in a uniform liquid state with a density $\rho_{0}=0.844$ and a temperature $T_{0}=0.722$, which is near the triple point of the LJ potential. This initial state of the LJ fluid is produced by a long-time (i.e., 2$\times 10^{7}$ time steps) quiescent MD simulation of the system shown in Fig. 1(a). The radial distribution function (RDF) of the initial liquid state is shown in Fig. 4(b).

At time $t$=0, the upper wall starts to move from left to right with a speed $V_{w}$, and the wall-driven shear flow is produced in the fluid domain between the walls. We note that the temperature of the fluid domain may significantly increase due to the viscous heating when the wall speed is sufficiently large.

III Results

We consider various channel widths $H$=168, 252, 336, 420, and 504, whereas the side lengths of the simulation box $L_{x}=L_{z}=16.8$ and the thickness of the wall $W\simeq 4.1$ are fixed. The speed of the upper wall also varies as $V_{w}=$1.0, 1.5, 1.75, 2.0, 2.5, and 3.0 for each channel width $H$.

The MD simulations are performed using the LAMMPS software package art:95P ; LAMMPS , in which Eq. (2) is time-integrated via the velocity Verlet method with a time-step size of $\Delta t=0.005$.

The width of the channel except the thin buffer layers along the boundaries (see Fig. 1(b)), i.e., $\frac{d}{2}<y<H-\frac{d}{2}$, is uniformly divided into 20 bins, and in each bin, the local macroscopic quantities, i.e., the density $\rho$, the flow velocity $u_{x}$, the temperature $T$, and the stress $p_{\alpha\beta}$, are calculated via the following equations.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the summation $\sum_{i}$ is taken over all the molecules and $V_{\rm{bin}}^{l}$ and $|V_{\rm{bin}}^{l}|$ represent the region of the $l$th bin and the volume of the local bin, respectively. On the right-hand side of Eq. (8), $\bm{f}_{i}$ is the force applied to the $i$th molecule due to the interaction among the ambient molecules, i.e., the right-hand side of Eq. (2).

The local quantities are also time-averaged in the stationary state after a long time $t_{0}$ has passed (i.e., $t_{0}=4\times 10^{7}\Delta t$), where the instantaneous quantities are sampled every 10 time steps in the interval $t=[t_{0},t_{0}+10^{5}\Delta t]$ (i.e., $10^{4}$ samples are averaged for each local quantity). The standard deviations of the instantaneous local macroscopic quantities shown in Fig. 2 are at most 0.023 for the velocity $u_{x}$, 0.018 for the temperature $T$, 0.0044 for the density $\rho$, and 0.15 for the normal stress $p_{yy}$.

In this section, we mainly show the results for $H$=504. The results for other channel widths are given in the supplemental materials.

III.1 Distribution of macroscopic quantities

Figure 2 shows the spatial distributions of the local macroscopic quantities (i.e., velocity, temperature, density, and stress) and local viscosity for different wall speeds $V_{w}$. The local viscosity is calculated as $\mu=p_{xy}/(du_{x}/dy)$. It is seen that the normal and shear stresses, $p_{yy}$ and $p_{xy}$, are uniformly distributed between the upper and lower walls for all cases. This fact confirms that the local stresses are balanced so that the flow velocity is in the stationary state.

The other macroscopic quantities spatially vary between the walls. The temperature increases in the middle region due to viscous heating, while it remains close to the wall temperature near the walls. By contrast, the local density decreases in the middle but increases near the walls.

Remarkably, for $V_{w}=$2.5 and 3, we can observe significant jumps in local density and local viscosity near the walls. Related to the rapid increase in the local viscosity in the vicinity of the wall, the velocity profile becomes nonlinear; i.e., the velocity gradient becomes much smaller near the wall than in the middle of the channel.

In the following text, we put focus on the peculiar behavior observed in the vicinity of the wall when the wall speed is large.

III.2 Solidification

Figure 3 shows snapshots of the local distributions of molecules in the lower, middle, and upper regions at two different time steps. From the motions of tracer particles, we can find the diffusive behaviors of local molecules in the different regions. It is clearly seen that the tracer molecules in the vicinity of the walls do not diffuse in the lateral direction ($y$-axis) but rather form a crystal-like structure, which is similar to that of the wall molecules. On the other hand, the molecules in the middle region diffuse in the lateral direction, as is observed in the fluid phase.

We can also observe the band structure in the solidified layer in Figure 3, where the bands run diagonally right upward near the upper wall. We note that the band structure can be observed when the wall speed and channel width are sufficiently large, e.g., $V_{w}\gtrsim 2$ and $H\gtrsim$ (see the supplemental materials). However, the direction of the band is not always diagonally right upward; instead, it may be the opposite direction or even appear in the $yz$–cross section. This finding indicates that band formation is not directly related to the flow velocity but rather to the compression of molecules from the bulk region toward the interface of the wall.

It is also seen from the supplemental materials that the solidification of the LJ fluid near the wall occurs only when the wall speed is sufficiently large, e.g., $V_{w}\geq 2$. When the wall speed is small, we observe only a thin absorption layer of molecules on the surface of the wall. The thickness of the absorption layer is only the length of a few molecules.

We also calculate the local RDF in the $l$th bin by

[TABLE]

where $n^{i}(r)$ counts the number of molecules within the distance [$r$,$r+dr$] from the $i$th molecule in the $l$th bin and $\left<\quad\right>_{l}$ represents the ensemble average over all the molecules contained in the local bin. Figure 4 shows the local RDFs of the molecules in the upper, lower, and middle regions. It is clearly seen that the RDFs near the walls have similar peak profiles to those of the FCC lattice structure of the wall, while the RDF in the middle region remains in the initial fluid state. This result also quantitatively confirms that the solidification occurs near the walls and that the lattice structures in the solidified layers are similar to those of the walls. In Section IV.2, we will discuss the effect of the wall structure in more detail.

To distinguish between the fluid and solidified phases, we measure the local lateral diffusion coefficient $D_{y}$ defined by

[TABLE]

where $v^{i}_{y}$ represents the lateral velocity of the $i$th molecules in the $l$th local bin and $\overline{\cdots}$ represents the time average in $\tau$.

Figure 5 shows the spatial distributions of the local lateral diffusion coefficient for different wall speeds. It is clearly seen that the lateral diffusion coefficient $D_{y}$ is negligibly small in the solidified or absorption layer near the wall. The solidified layer with a very small lateral diffusion coefficient, e.g., $D_{y}<0.01$, broadens when the wall speed changes from $V_{w}=$1.75 to 2.0.

We summarize the results of the solidification under different parameters in Fig. 6, in which a diagram of the solidification vs. the channel width $H$ and the square of the wall speed $V^{2}_{w}$ is shown. In the figure, the thickness of the solidified layer $l_{S}$, which is defined by the thickness of the layer where the local lateral diffusion coefficient is smaller than 0.01, i.e., $D_{y}<0.01$, is indicated by the color legend. The square symbols $\square$ represent the results obtained when the solidified layer extends far beyond the molecular size, i.e., $l_{S}>10$. It is clearly seen that broadened solidification occurs only when the wall speed is sufficiently large, e.g., $V_{w}\geq 2$. In the high-speed regime, $V_{w}\geq 2$, the thickness of the solid layer $l_{S}$ proportionally increases with the channel width $H$ (see Figure 7). The mechanism underlying these observations will be discussed in the next section.

Incidentally, in Figure 6, the square of the wall speed $V_{w}^{2}$ is used as the horizontal axis rather than the wall speed itself. This is because $V_{w}^{2}$ represents the amplitude of the viscous heating relative to the thermal conduction in macroscopic energy transport, i.e.,

[TABLE]

where $\lambda$ is the thermal conductivity and $\Delta T$ is a characteristic temperature rise.

IV Discussion

IV.1 Hydrodynamic Explanation

In this subsection, we consider the mechanism of the solidification from a hydrodynamic point of view. We suppose a Newtonian fluid with a constant viscosity $\mu$ and the Fourier law of heat conduction with a constant thermal conductivity $\lambda$. We also introduce the normalized coordinate $\hat{y}=y/H$, which is relevant to the hydrodynamic analysis. In the following part of this section, we consider only the stationary state. Then, the spatial distribution of temperature is described by

[TABLE]

with the boundary condition $T=T_{w}$ at $\hat{y}=$0 and 1. The solution to the above equation is explicitly calculated as

[TABLE]

where $\Delta T=\mu V_{w}^{2}/8\lambda$ is the difference in temperatures between the region at the wall and the region in the middle of the channel.

The mass conservation is written as

[TABLE]

where $\rho_{0}$ is the initial density of the LJ fluid. We can also easily obtain from the momentum balance equation that the bulk pressure is spatially uniform in the stationary state because of the continuity condition $\frac{\partial v_{y}}{\partial\hat{y}}$=0.

We suppose that the equation of state $\rho={\cal F}(p,T)$ holds at the local fluid elements in the stationary state even under shear flow and satisfies the conditions $\left(\frac{\partial{\cal F}}{\partial T}\right)_{p}<0$ and $\left(\frac{\partial{\cal F}}{\partial p}\right)_{T}>0$. Then, the bulk pressure $p$ is determined from the equation

[TABLE]

This indicates that the bulk pressure does not depend on the channel width $H$ but depends only on the wall speed $V_{w}$ when the initial states $\rho_{0}$ and $T_{w}$ are fixed. Note that $T(\hat{y})$ does not depend on the channel width $H$ in Eq. (13).

The dependency of the bulk pressure on the wall speed is obtained by taking the derivative of Equation (15) against $V_{w}$, i.e.,

[TABLE]

Note that from Equation (13), $\frac{\partial T}{\partial V_{w}}>0$ holds at any local position $\hat{y}\in(0,1)$. Thus, the bulk pressure dose not depend on the channel width $H$ but monotonically increases with the wall velocity $V_{w}$. In fact, in Figure 8, our simulation results demonstrate that the bulk pressure monotonically increases with the wall speed but is less dependent on the channel width.

The local density is described by the equation of state as $\rho(\hat{y})={\cal F}(p,T(\hat{y}))$, where the bulk pressure $p$ is spatially uniform and the temperature is described by Equation (13). It is seen that the local density monotonically increases while approaching the wall (under the condition $\left(\frac{\partial{\cal F}}{\partial T}\right)_{p}<0$) and takes the maximum value at the wall, i.e., $\rho(\hat{y})\nearrow\rho_{w}={\cal F}(p,T_{w})$ as $\hat{y}\rightarrow$ 0 or 1. This result indicates that the maximum of the local density also monotonically increases with the wall speed $V_{w}$ but does not depend on the channel width $H$ in the same way as the bulk pressure.

If we suppose that the fluid is solidified when the local density exceeds a critical density $\rho^{*}$, then the dependency of the maximum density on the wall speed and channel width indicates that the solidification never occurs unless the wall velocity exceeds the critical velocity $V_{w}^{*}$, which is obtained from $\rho^{*}={\cal F}(p(V_{w}^{*}),T_{w})$, regardless of the channel width $H$. This concisely explains the observation of the existence of a critical wall speed for solidification in Figure 6.

The reason why the thickness of the solidified layer, $l_{S}$, is proportional to the channel width $H$ is also explained; i.e., the local density is a monotonic function of $\hat{y}$ and independent of $H$, so the solidified layer, where the local density is larger than the critical density, i.e., $\rho(\hat{y})>\rho^{*}$, is uniquely determined by the condition $\rho(l_{S}/H)=\rho^{*}$ for a given wall velocity. Thus, the thickness of the solidified layer $l_{S}$ is proportional to the channel width $H$.

Figure 9 shows the transient behaviors of the local $\rho$–$T$ states between the bulk ($y=60$) and the interface of the wall ($y=0$). Here, instead of using Eq. (5), we calculate the local density by

[TABLE]

where the local RDF $g^{l}(r)$ is defined in Eq. (9). It is seen that the conditions supposed above, i.e., $(\frac{\partial{\cal F}}{\partial T})_{p}<0$ and $(\frac{\partial{\cal F}}{\partial p})_{T}>0$, are relevant to the simulation results and, in fact, the maximum density increases with the wall speed $V_{w}$.

The local density rapidly increases around the solidification line while approaching to the wall from the bulk, and the $\rho$–$T$ state enters into the liquid/solid coexistence regime in the phase diagram. For a large wall speed $V_{w}=2.5$, remarkable solidification (or crystallization) is observed even in the liquid/solid coexistence regime, where the thermal expansion becomes very small, i.e, $|\frac{1}{\rho}(\frac{\partial\rho}{\partial T})_{p}|\ll 1$, as is usually observed in solid materials. For a small wall speed $V_{w}<2$, solidification is remarkably not observed, but a thin absorption layer forms on the surface of the wall (see also Fig. 5).

These observations seem to indicate that the LJ molecules in the fluid phase are confined in the vicinity of the wall due to the thermohydrodynamic coupling and that when the density in the confined regime is close to the solidification line, the tightly confined LJ molecules are solidified via the interaction with the wall molecules. This also indicates that both the wall structure and the $\rho$–$T$ state in the vicinity of the wall strongly affect the solidification near the wall.

IV.2 Effects of the wall structure and the initial state

Thus far, we have considered the face-centered-cubic (FCC) lattice structure for the wall molecules and the initial condition of the LJ fluid near the triple point, i.e., $\rho_{0}$=0.844 and $T_{0}$=0.722. In this subsection, we change the wall structure and the initial condition and investigate the effects of the wall structure and the initial condition of the LJ fluid.

Figure 10 shows snapshots of the molecules composing the bottom wall (i.e., $-H<y<0$) and those in the vicinity of the wall (i.e., $0<y\leq 60$) for the channel width $H\simeq 500$ and the wall velocity $V_{w}$=3.0. It is seen that solidification occurs with both the FCC and BCC structures.

However, interestingly, the lattice structures of the solidified layer are different from each other. With the BCC wall, the LJ molecules in the solidified layer also create the BCC lattice structure, although it is known that the FCC structure appears during the crystallization of the LJ molecules in the equilibrium state. The thickness of the solidified layer is also affected to the wall structure; i.e., the solidified layer for the BCC wall is thinner than that for the FCC wall.

Figure 11 shows the transient behaviors of local $\rho$–$T$ states from the bulk regime to the surface of the wall for three different initial states of the LJ fluid. If the initial state is close to the solidification line (i.e., the square $\Box$), remarkable solidification is observed in the liquid/solid coexistence regime in the phase diagram. Even if the initial state is slightly away from the solidification line (i.e., the circle $\bigcirc$), we can observe that the solidified layer forms in the vicinity of the wall. However, when the initial state is far from the solidification line (i.e., the triangle $\bigtriangledown$), the local $\rho$–$T$ state cannot approach the solidification line even in the vicinity of the wall, so solidification does not occur in the vicinity of the wall.

From these observations, we can conclude that the LJ molecules are confined in the vicinity of the wall via thermohydrodynamic transport and that when the local $\rho$–$T$ state is close to the solidification line in the vicinity of the wall, the LJ molecules are solidified due to the interaction with the crystallized wall molecules.

V Concluding remarks and perspectives

We carried out MD simulations of the thermohydrodynamic lubrication of the LJ fluid between atomistic thermal walls. A counterintuitive solidification, i.e., viscous heating-induced solidification, was discovered, in which the LJ fluid is solidified near the wall because of the viscous heating generated in the bulk regime. It was found that the solidification occurs only when the wall speed is sufficiently large regardless of the channel width, even though the thickness of the solidified layer increases with the channel width. Band formation was also found in the solidified layer when the channel width is large.

We investigated the solidification mechanism in detail from both macroscopic and microscopic points of views. It was found that the LJ molecules are confined in the vicinity of the wall via thermohydrodynamic transport and that when the local $\rho$–$T$ state in the vicinity of the wall is close to the solidification line in the phase diagram, the LJ molecules are solidified due to the interaction with the crystallized wall molecules. The lattice structure of the wall molecules affects that in the solidified layer near the wall; thus, the BCC lattice structure is created in the solidified layer even though the FCC lattice structure is more stable for the LJ molecules.

This study explicitly demonstrates that even for a simple fluid in a simple geometry with smooth boundaries, the thermohydrodynamic coupling in high-speed lubrication flow and the molecular interaction at the interface produce unexpected flow behavior. In the literature art:13Aetal , it is reported that nanoscale surface texture of wall and molecular structure of fluid significantly affect the solidification in confined geometries. Investigation on the effects of the roughness of wall surface and the molecular structure of fluid should represent an important future research direction.

This study also illustrates that the molecular interaction between wall and fluid significantly affects the phase transition behavior in the confined regime near the wall. This result gives an important message for the future development of multiscale simulation technologies for large-scale complex flows.

Acknowledgements

This study was financially supported by JSPS KAKENHI Grant Number 16K17554 and 17H01083.