Finite-Size Effects on Karman Vortex in Molecular Dynamics Simulation

TL;DR

This paper investigates how finite-size effects influence the simulation of Karman vortices in molecular dynamics, aiming to identify conditions for quantitatively matching experimental vortex shedding frequencies.

Contribution

It identifies key factors like vortex interference and viscosity nonuniformity that cause finite-size effects in MD simulations of Karman vortices.

Findings

Finite-size effects significantly alter vortex characteristics.

Simulation conditions can be optimized for quantitative accuracy.

Finite-size effects are mainly due to vortex interference and viscosity nonuniformity.

Abstract

The characteristics of the Karman vortex generated by a molecular dynamics (MD) simulation exhibit strong finite-size effects, and MD can only reproduce the experimental results qualitatively. Here, we seek the simulation conditions for quantitatively reproduce the vortex shedding frequency. We found that the finite-size effects are mainly caused by the interference of vortices and the nonuniformity of the fluid viscosity coefficient.

Equations

$$St=0.212-\frac{4.5}{Re},$$

Full text

It is well known that an obstacle in a flow generates two rows of staggered vortices called Kármán vortex street behind the obstacle. Flows around an obstacle can be found in many places around us, such as airflow around a building, water flow around a screw propeller, and so forth. Because the Kármán vortex is a source of noise and vibration, it is crucially important to understand the flow characteristics around an obstacle in the engineering field. Especially, the flow around a circular cylinder is one of the most fundamental problems of fluid dynamics. Therefore, experimental and numerical studies have been conducted on the flow around the circular cylinder williamson96 . The shedding characteristics of the vortices are phenomenologically well understood for the Newtonian fluid. The shedding frequency is given experimentally by the following equation roshko54 :

[TABLE]

where $St=fD/V$ and $Re=\rho DV/\eta$ are the Strouhal number and Reynolds number, respectively ($f$: vortex shedding frequency, $D$: cylinder diameter, $V$: inflow velocity, $\rho$: fluid density, $\eta$: fluid viscosity coefficient).

With the progress in computational power, Kármán vortex can be generated by molecular dynamics (MD) simulation rc86 ; rapaport87 ; awn18 . Although MD simulation is advantageous to deal with the complex flow phenomena directly, such as the Toms effect gadd65 exhibited by polymer solutions, and cavitation in multiphase flow gm16 , it suffers from strong finite-size effects. To understand the finite-size effects on MD, we investigate the difference between the Strouhal number $St^{\rm(MD)}$ obtained by MD simulation and the value predicted by Eq. (1). This relation Eq. (1) holds for not only macroscopic methods based on the Navier-Stokes equations but also mesoscale methods hd97a ; lg02 . However, the values obtained by the MD simulations are twice as large as the expected value rc86 ; rapaport87 ; awn18 . Because both $St$ and $Re$ are dimensionless quantities that do not depend on a typical spatial scale such as the cylinder diameter, it is expected that the relation Eq. (1) also holds in MD simulation. Possible causes of the deviation are the influence of compressibility sd03 ; ouv06 due to high flow velocity and the interference of vortices due to periodic boundary conditions swp99 ; twt14 . The effects can be investigated by changing the density for the compressibility and the channel width for the vortex interference. Therefore, in the present study, we investigate the condition that quantitatively reproduces Eq. (1) for the Kármán vortex behind the cylinder by the MD simulation.

We used the Weeks–Chandler–Andersen potential wca71 $\phi=4\epsilon\left[\left(\sigma/r\right)^{12}-\left(\sigma/r\right)^{6}+1/4\right]\theta(2^{\frac{1}{6}}-r)$ for the interparticle interaction of the fluid, where $\theta$ is the Heaviside function, $r$ is interparticle distance, $\epsilon$ and $\sigma$ represent the energy and the length scale, respectively. The mass of the particle is $m$. Hereafter, physical quantities are expressed in units of energy $\epsilon$, length $\sigma$, and time $\tau_{0}=\sigma\sqrt{m/\epsilon}$. The simulation box is almost the same as in Ref. awn18 and is a rectangle with dimension $L_{x}\times L_{y}$, where $L_{x}=3600$ and $500\leq L_{y}\leq 5000$. The flow direction of the fluid is the $x$-direction. The periodic boundary condition is taken for all directions. The circular cylinder is modeled by a set of particles whose positions are fixed on the cylinder surface. The cylinder with diameter $D=100$ is located at $(x,y)=(1000,L_{y}/2)$. To control the temperature and the inflow velocity, we used the Langevin thermostat awn18 in the region of $3200\leq x\leq 3600$. The friction coefficient of the Langevin thermostat ($\zeta$ of Eq. (4) in Ref. awn18 ) is $0.1$ and is linearly increased from $0.0001$ to $0.1$ at $3200\leq x\leq 3400$. The temperature is set to $T=1$. This thermostat relaxes the velocities of the fluid particles into the Maxwell distribution whose average velocity is given by $V$ in the $x$-direction. The density is in the range of $0.4\leq\rho\leq 0.83$. In the three-dimensional (3D) simulations, the thickness of the simulation box in the $z$-direction $L_{z}=10$ is employed. When Kármán vortices appear, a periodic lifting force $F_{\rm L}$ whose period is equal that of the vortex shedding acts on the cylinder. Here, we adopted the frequency of lift coefficient $C_{\rm L}=2F_{\rm L}/(\rho DV^{2})$ as the shedding frequency of the Kármán vortex. MD simulations were performed using the velocity-Verlet algorithm with a time step of $0.004$, for up to a maximum of $10~{}000~{}000$ time steps. The time integration is performed by LAMMPS plimpton95 .

Figure 1(a) shows the $Re$ dependence of $St^{\rm(MD)}$ in the two-dimensional (2D) simulations. The broken line in the figure shows Eq. (1). While all results deviate from Eq. (1) in the case of $L_{y}=500$, the results of $L_{y}=5000$ tend to approach Eq. (1). Because we controlled $Re$ by changing $V$, the Mach number $Ma$ increases with increasing $Re$. The decreasing in $St$ of $L_{y}=5000$ at high Reynolds number ($Re\gtrsim 160$) is due to the high Mach number $Ma$. Because $Ma$ is 0.5 or more in the range of $Re\gtrsim 160$, the compressibility affects the vortex shedding frequency. The compressibility effects on the frequency also appear in the simulation based on the Navier-Stokes equations sd03 ; ouv06 . In addition, because the flow velocity is too fast in the MD simulation, the formation of vortices is inhibited due to the detachment of the fluid behind the cylinder. The inset of Fig. 1 shows the $L_{y}$ dependence of the relative error of $\rho=0.40,0.83$ at $Re=100$. The relative error is defined as $\left|St^{\rm(MD)}-St\right|/St$. In the case of $\rho=0.40$, the relative error decreases monotonically as $L_{y}$ increases. Therefore, the interference of the vortex with the image cell is one of the main causes of the deviation from Eq. (1). As for $\rho=0.83$, the relative error only reaches up to $\sim 30\%$ and the relative errors are almost unchanged between $L_{y}=2500$ and $L_{y}=5000$. This observation suggests that another cause also exists for this phenomenon at the high density. Figure 1(b) shows the results of the 3D simulations. Like the 2D simulations, $St^{\rm(MD)}$ tends to approach Eq. (1) as $L_{y}$ increases.

The density dependence of the relative error at $Re=100$ is shown in Fig. 2(a). As the density decreases, the error at $L_{y}=5000$ decreases monotonically while it is almost independent at $L_{y}=500$ for the low density ($\rho\leq 0.6$). Therefore, the deviation from Eq. (1) decreases as the density decreases when the system is almost free from the interference of the vortices. As shown in Fig. 2(b), the viscosity increases rapidly with the density. Therefore, the variation in the viscosity due to the density change also increases with increasing the density. In flows, the spatial and temporal density variation is caused by the pressure change. Figure 2(c) shows the time-averaged local density $\rho_{\rm l}$ and the temporal density fluctuation $\Delta\rho$ downstream of the cylinder at Re = 100 with $L_{y}=5000$ in the 2D case. A characteristic change in $\Delta\rho$ exists at $200<\Delta x<300$ in the case of $\rho=0.83$. The maximum density fluctuation of $\rho=0.40$ and $\rho=0.83$ are $\Delta\rho=0.035$ and $\Delta\rho=0.12$, respectively, where $\Delta x$ is the downstream distance from the center of the cylinder. Since $Ma$ is not so different ($Ma\simeq 0.4$ and $0.5$ at $\rho=0.40$ and $0.83$, respectively), it is not caused by $Ma$ change. The average density behind the cylinder is spatially lower than the other regions at higher density. Since these density variations induce nonuniform viscosity, the effective Reynolds number which determines the vortex shedding frequency is altered. Hence, we conclude that the nonuniformity of viscosity near the cylinder is the cause of $St$ deviation at high density. The relative errors in 3D are smaller than those in 2D because the viscosity variation is smaller as shown in Fig. 2.

In summary, we evaluate the Strouhal number of a Kármán vortex behind a circular cylinder by MD simulations. The main causes of the deviation from the experimental results are the interference of vortices and the nonuniformity of the viscosity. The former can be suppressed by widening the channel width, and the latter can be reduced by decreasing density. Therefore, we conclude that a nanoscale Kármán vortex can be analyzed quantitatively by the MD simulation.