Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Large-scale molecular dynamics simulations of Kármán
vortex and sound wave: cavitation and polymer effects
Yuta Asano 1, Hiroshi Watanabe 2, and Hiroshi Noguchi 1
1

Institute for Solid State Physics,

The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581, Japan2

Department of Applied Physics and Physico-Informatics,

Keio University, Yokohama, Kanagawa 223-8522, Japan

Abstract

lives. Therefore, a thorough understanding of the

In many industrial flows, microscale effects

flow phenomena will improve our lives, also

have a significant impact on their macroscalesignificantly

flows. Turbulent drag reduction by polymer

environmental protection. Flow control in

addition and cavitation around rotating turbines

industrial processes is an essential topic in fluidin water are typical examples. In these flows,engineering, such as turbulence drag reduction

phenomena in a wide range of scales coexist,

by polymer addition [5] and gas-liquid multi-

from molecular scale motion to macroscale fluid

phase flow due to the cavitation [6]. Thus, many

motion, and these interact with each other underresearchers have conducted experiments andstrongly non-equilibrium. Therefore, it isnumerical simulations based on the Navier—challenging to clarify the underlying mechanism

Stokes equation to clarify the mechanisms.

by its complexity.

However, since various phenomena in a wide

We have performed large-scale molecularrange of scales coexist and interact with each

dynamics (MD) simulations of complex flows,other in non-equilibrium, the analyses are

such as polymer solutions and cavitating flow, as

challenging issues. Therefore, it is crucial towell as MD of sound propagation, using theanalyze all phenomena simultaneously, from

ISSP supercomputer. In this activity report, wemolecular-scale dynamics to fluid motion topresent a brief overview of these studies [1-4].

understand the complex flows.

contributes

toglobal

Mesoscale analysis methods such as the

1. Introduction

lattice Boltzmann method [7] and Multi-particle

Fluid flows play a vital role in a wide range

collision dynamics [8,9] have been developed toof fields, from tap water familiar to our daily

analyze the flow considering the effects oflives to fuel pipelines in energy supply systems.

molecular-scale dynamics. Also, the micro-

It is no exaggeration to say that flow sustains our36Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

macro coupling approaches [10,11] have been
proposedtoinvestigate

more

detailed

molecular-scale information. However, in lightof today's computational power, it is possible tosimulate flows by all-particle simulations,
namely a molecular dynamics (MD) simulation.
In the MD simulation, molecular-scale dynamics,
such as phase transition and polymer motion,
naturally appear in the flow by simply solving
the Newtonian equations of motion, only
providing the interparticle interactions. In
conventional

fluid

analysis

methods,

phenomenological models, such as constitutive
equations and an equation of state, are required
to describe these phenomena. Instead, we haveto prepare a sufficiently large computational cell
that

can

resolvecharacteristics.

themacroscale

flow

Therefore,

althoughitisinevitably a large-scale simulation, it will help

Fig. 1: Instantaneous vorticity fields for

us understand complex flow phenomena.

(a) the Newtonian fluid, (b) polymer

In this report, we present a brief overview offlow

analysisbyusing large-scale

solution with the number of segmentsMD

𝑵𝐬 = 𝟏𝟎, and (c) polymer solution withsimulations. We explain the effect of polymer

𝑵𝐬 = 𝟏𝟎𝟎 at the Reynolds number 𝑹𝒆 =addition and cavitation on the flow around a𝟔𝟒 [1].

cylinder in sections 2 and 3, respectively. In
Adding a tiny amount of polymer into a

section 4, we describe the application to

Newtonian fluid involves significant changes insoundwave propagation in a simple fluid. In

the flow characteristics, such as a dramatic

section 5, we explain the recent study of

reduction in frictional drag [5] and suppression

cavitating flow in a polymer solution. Section 6of the Kármán vortex [12]. Since most industrial

is a summary.

flows are turbulent, the flow control by adding
polymers is an attractive method in engineering

2. Effects of polymer addition on

applications

Kármán Vortex [1]37

becauseitsaves

energyActivity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

necessary to elucidate the effects of polymers on
vortices from the molecular scale. We focus on
the flow around a cylinder, which is the mosttypical flow accompanied by vortices and issuitable for analysis due to its periodicity. We
investigate the effects of polymers on vortices inthe two-dimensional flow around a circularcylinder by MD simulations.
Figure 1 shows typical snapshots of thevorticity field. As for the Newtonian fluid, analternating vortex street, the Kármán vortex
street, appears behind the cylinder. Similar
vortices appear for the short-chain polymers.
However, for long-chain polymers, the vortices
are blurred. Therefore, the chain length of thepolymer significantly affects vortex formation.
Fig. 2: Fourier spectra of the lift

Figure 2 shows the Fourier spectra of the lift

coefficient at 𝑹𝒆 = 𝟔𝟒 [1]. (a) Newtonian

coefficient. For the Newtonian fluid, a sharp

fluid. (b) Polymer solution with 𝑵𝐬 = 𝟏𝟎

peak appears, which indicates periodic vortex

and volume fraction 𝝓 = 𝟏𝟎. (c) Polymer

shedding. Note that the peak frequency issolution𝝓=slightly increased by the finite size effects we

𝟎. 𝟎𝟐.(d) Polymer solution with 𝑵𝐬 = 𝟏𝟎𝟎

reported in Ref. 13. The short-polymer solutionand 𝝓 = 𝟎. 𝟏𝟎𝟕 .

andshows a similar sharp peak. On the other hand,horizontal dotted lines denote the peak

the peak position shifts to a lower frequency, andposition and height for the Newtonian

the spectrum broadens for the long-polymer

fluid, respectively.

solution. These effects enhance with increasingwith

𝑵𝐬 = 𝟏𝟎𝟎
Theandverticalpolymer concentration. The effects of polymer
consumption. However, the mechanism is not

addition on the flow, (1) vortex blurring, (2)well understood because the polymer rheology

frequency reduction, and (3) chain length effect,in a flow is unclear.

obtained by MD simulations, are in good

The interaction between polymers and vortices

agreement with experimental results [12].

presumably plays a significant role in the

Figures 3(a) and 3(b) show the gyration

changes in flow characteristics. Therefore, it isradius of short- and long-chain polymers,

38

Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

the other hand, long-chain polymers are
significantly stretched near the cylinder and in
the vortex behind the cylinder. Furthermore, thevortex also entrains the polymer in its vicinity.
Figures 3(c) and 3(d) show the orientational
order Q of the short- and long-chain polymers,
respectively. When the polymer aligns in theflow direction, 𝑄 = 1

the polymer aligns

perpendicular to the flow, 𝑄 = −1 the polymer
directs randomly, 𝑄 = 0 . For the short-chain
polymer, polymers orient randomly in the entire
area. However, for long-chain polymers, thepolymers are aligned in the flow direction near
the cylinder. Furthermore, the polymers also
tend to align in the flow direction in the vortices
behind the cylinder. The reduction of theorientational order behind the cylinder is due tothe entrainment by the vortex. Therefore, thestretching and entrainment of the polymer play asignificant role in vortex suppression. As aconsequence, the shapes of the vortices become
blurred.
Fig. 3 Spatial distributions of the gyration

3. Effects of cavitation on the Kármán

radius 𝑹𝐠 and orientational order Q for
the polymer solutions with (a, c) 𝑵𝐬 = 𝟏𝟎
and (b, d) 𝑵𝐬 = 𝟏𝟎𝟎 at 𝑹𝒆 = 𝟔𝟒 and 𝝓 =𝟎. 𝟏𝟎𝟕[1]. 𝑹𝐠𝟎 is the gyration radius in theabsence of the flow.

vortex behind circular-cylinder arrays[2]
Cavitation is a phenomenon of forming
bubbles due to local pressure drop in a flowing
liquid with high velocity [14]. When cavitation

respectively. For short-chain polymer, theoccurs in fluid machinery, such as pumps andgyration radius is virtually identical to that in theturbines, it causes various adverse effects suchstationary fluid. Therefore, the polymer flows

as performance degradation, noise and vibration,with an almost spherical shape as a particle. On

and erosion. Therefore, it is vital in fluid

39

Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Fig. 5: Instantaneous void fraction fieldsfor (a) 𝑻 = 𝟏. 𝟑 , (b) 𝑻 = 𝟏. 𝟐𝟓 , and (c) 𝑻 =𝟏. 𝟐 [𝟐].

Fig 4: Instantaneous vorticity fields forthe (a) non-cavitating flow (𝑻 = 𝟐) andthe vorticity and void fraction field, respectively.

(b-d) cavitating flows (𝑻 = 𝟏. 𝟓, 𝟏. 𝟐𝟓,

At high temperature (𝑇 = 2), no bubbles appear

and 1.2) [2].

(i.e., non-cavitating flow). The Kármán vortex
appears behind the cylinders, and the upper andengineering to clarify the cavitation mechanism

lower vortex streets synchronize in anti-phase.

for appropriate control of the cavitating flow.
However,

becauseofthelimitations

As temperature decreases, bubbles appear in theofvicinity

discussing the dynamics of microscopic bubble

of

the

cylinderat 𝑇 = 1.3 inconjunction with vortex shedding (Fig. 5 (a)).

nuclei in a macroscopic flow, the cavitationThe

mechanism is still not understood well. To

vortexstreets

behindthecylindersynchronize in anti-phase, as in the case of 𝑇 =clarify the effects of cavitation on the flow, we

2 (Figs. 4(a) and 4(b)).

performed a large-scale MD simulation of the

At 𝑇 = 1.25, the gas phase region is formed

flow for the Lennard-Jones (LJ) aroundbehind the cylinder and is attached to thecylinders in a side-by-side arrangement.

cylinder (Fig. 5(b)). Also, the phase difference

Figures 4 and 5 show the typical snapshots of

between the upper and lower vortex streets40

Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyoof the lift coefficient 𝐴L and the vortex
formation length 𝐿f . For comparison, the figure
also shows the results for the Newtonian fluid.
For the reference liquid, 𝐴L and 𝐿f tend toincrease and decrease with increase Re,
respectively. However, in the case of LJ fluids,
they exhibit non-monotonic behavior to 𝑅𝑒 .
This behavior is probably due to the conflicting
effects of increasing 𝑅𝑒 and bubble generation.
𝐴L rapidly decreases and finally disappears dueto the bubble generation. Simultaneously withthe disappearance of 𝐴L, 𝐿f rapidly increases.
For the Newtonian fluid, 𝐴L and 𝐿f show

Fig. 6: 𝑹𝒆 dependence of (a) the lift

almost the same trend for the two-cylinder

coefficient amplitude and (b) vortex

system and the one-cylinder system, which hasformation for the LJ fluid and Newtonian

only one cylinder with the same distance

fluid [2]. Double and Single in the legend

between the neighboring (periodic image)denote two- and one-cylinder systems,

cylinders (therespectively.

vortex streets are always

synchronized in-phase). In the case of LJ fluid,slightly shifts (Fig. 4(c)). At 𝑇 = 1.2 , the gas

however, the behavior before bubble generation

phase region behind the cylinder expands further

is different. The upward and downward peaks(Fig. 5. (c). In addition, the upper and lower

appear in 𝐴L and 𝐿f , respectively. Theirvortex structures become asymmetric (Fig. 4(d)).

positions for one- and two-cylinder systems areThis asymmetric vortex structure switches over

𝑅𝑒 = 105 and 101, respectively. These peaks

a long period. For the fluid without phaseare probably due to a small nucleus generationtransition, such a change in the flow field does

which is the precursor of the cavitation or

not occur with temperature change. We founddensity fluctuations originating from the criticalthat the density fluctuation caused by the bubble

phenomena. Since the phase retraction amplifies

generation has a significant effect on the vortex

the oscillation amplitude, the peak appears at astructure.

higher temperature in the two-cylinder system

As the vortex structure changes, the lift force

than in the one-cylinder system. After bubble

acting on the cylinder is also affected. Figure

generation, the bubbles weaken the interference6(a) shows the 𝑅𝑒dependence of the amplitude

between the vortex street, and the effect of the41

Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Fig. 8: Waveforms of the LJ fluid for theFig 7: Waveforms of the LJ fluid for the(a) 𝒇 = 𝟎. 𝟎𝟎𝟐 and 𝑨 = 𝟓, (b) 𝒇 = 𝟎. 𝟎𝟎𝟐(a) amplitude 𝑨 = 𝟐. 𝟓 (b) 𝑨 = 𝟓 (c) 𝑨 =and 𝑨 = 𝟏𝟎 (c) 𝒇 = 𝟎. 𝟎𝟎𝟒 and 𝑨 = 𝟓 ,𝟏𝟎 , and (d) 𝑨 = 𝟐𝟎 at frequency 𝒇 =and (d) 𝒇 = 𝟎. 𝟎𝟎𝟒 and 𝑨 = 𝟏𝟎 [𝟑] . The

𝟎. 𝟎𝟎𝟏 [3]. The red and black lines

red and black lines represent the results

represent the results of MD simulationof MD simulation and the Burgers'

and the Burgers' equation, respectively.

equation, respectively.
phase difference is almost negligible. We foundthat the bubble generation significantly affects

extend in a wide range of applications, such as

the vortex structures, vortex interference, and lift

medical treatment and food processing. For

due to vortex shedding.

these applications, it is crucial to understand thecharacteristics of ultrasonic cavitation and the

4. Soundwave propagation in simple

propagation in complex fluids. Therefore, it is

fluid [3]necessary

toanalyze

the

sound

wave

propagation from the molecular scale.

Sound waves are a familiar phenomenon in
our daily lives, such as the voices and sound of

We apply the MD simulation to the soundmusical instruments. Because of the soundpropagation in the LJ fluid. Moreover, weproperties of propagating a wide range of solid,

perform the fluid dynamics calculation utilizing

liquid, and gaseous materials, their application

Burgers' equation to justify the validity of theexists in a wide range of engineering fields. In

MD simulation for the sound waves. When theparticular, applications of ultrasonic cavitation

same fluid is analyzed, the two results should

42

Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyodeviatethe magnitude of the deviations

becomes greater for higher amplitudes. Thisdeviation is due to the acoustic flow that occurs
by the high nonlinearity owing to high frequency.
Since Burgers' equation is not applicable in theregion where acoustic flow occurs, a fluid
analysis

involvinghigher-order

terms

is

necessary for a quantitative discussion. The MD
simulation provides a waveform reflecting thenonlinearity, thus obtaining a
Fig.

9:

𝒙

dependences

of

the

frequency-

dependent sound speed.

(a)normalized amplitude and (b) attenuation

The attenuation coefficient of a sound wave

coefficient of the LJ fluid [3]. The solid

is estimated from the decrement in amplitude.

and broken lines represent the results of

Figure 9(a) shows the amplitude decrement at

MD simulation and Burgers' equation,

𝑓 = 0.001 . When the amplitude is small, therespectively. The dotted line in panel (b)deviation of the MD simulation from theshows

Burgers’ equation is slightly different due to thetheclassical

attenuation

thermal fluctuations. As the amplitude increases,

coefficient.

the agreement between them becomes better. We
agree within the statistical error.

estimate the attenuation coefficient by the

Figure 7 shows the amplitude dependence oflogarithmic decrement rate of the amplitude at

the waveform of the LJ fluid at frequency 𝑓 =each position 𝑥 . Figure 9(b) shows the0.001 . For small amplitudes, the waveform isattenuation coefficient as a function of 𝑥 . The

sinusoidal owing to the weak nonlinearity. As

dotted line shows the value of the attenuation

the amplitude increases, the waveform becomescoefficient evaluated from the classical theory.

a sawtooth waveform (Figs. 7(c) and 7(d)). The

The results of the Burgers’ equation approach

black dotted lines in Fig. 7 are the numericalasymptotically to the classical attenuation

solutions of Burgers' equation, which is in goodcoefficient at a sufficient distance from theagreement with the results of the MD simulation.

sound source position ( 𝑥 = 0 ). Although the

Figure 8 shows the frequency dependence of

results of the MD simulation are in good

the waveform at higher frequencies. As theagreement with the results of Burgers' equation

frequency increases, the waveforms of the MD

when the amplitude is large, the attenuation issimulation and Burgers' equation become

overestimated by the nonlinearity. On the other43

Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Fig. 10: Instantaneous void fraction fields of the (a) LJ fluid and (b) polymer solution, andtime-averaged void fraction of the (c) LJ fluid and (d) polymer solution at 𝑻 = 𝟏. 𝟐𝟓 [𝟒].
hand, the attenuation coefficient is hardly

flows, it has attracted attention. Although it isestimated for small amplitudes due to thermal

challenging to simultaneously analyze polymer

fluctuation. Therefore, to obtain the classicalmotion,attenuation coefficient by the MD simulation,

macroscopic flow, the analysis should be

the practical method is to estimate thepossible by integrating the methods described inparameters of the Burgers’ equation from theSections 2 and 3. Here we briefly describe only

waveform of MD simulation.

the main results.

bubblenucleusdynamics,

andWe investigate the effect of the polymer

5. Effects of the polymer addition on the

addition on the cavitating flow around a circular

cavitating flow [4]cylinder. Figure 10 shows the distribution of

Finally, we briefly describe our recent

void fraction in the cavitating flow of the LJ

analysis on the cavitating flow. In extreme cases,

fluid at 𝑇 = 1.25 , at which the cavitation occurs

the adverse effects on fluid machinery caused byas described in Sec. 3. The polymer addition

cavitation are prevented by controlling the flow

dramatically suppresses theto suppress cavitation. However, in general,

bubbles. Note that the polymer addition has

most industrial flows are highly irregularnegligible effects on the phase transition in the

turbulent fields, and it is difficult to achieve suchabsence of the flow. From the temperature

flows by mechanical design. Since polymer

changeaddition is one of the methods to realize suchconformation change in the polymer, we found

44

in

the

fluid

formation ofelement

andtheActivity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

that the vortex suppression and the entropic

Morii, Shinichi Tsuda, Yuichi Kunishima, andelasticity by the polymer have essential roles in

Yuji Higuchi for helpful discussions. This

the suppression of cavitation.

researchwassupported

by

MEXT

as

“Exploratory Challenge on Post-K computer”

6. Summary(Challenge

of

Basic

Science

Exploring

We performed large-scale MD simulations

Extremes through Multi-Physics and Multi-

for complex fluids involving polymers, gas-

Scale Simulations) and JSPS KAKENHI (Grant

liquid multi-phase flows owing to cavitation,No.

and sound waves. The direct analysis from theSupercomputer Center, Institute for Solid State

molecular scale enables us to analyze the flow

Physics (ISSP), University of Tokyo, Research

field reflecting the micro-scale effects, such asCenter for Computational Science (RCCS),

the elongation properties of polymers, phase

Okazaki, Japan, and Center for Computationaltransitions that are difficult to model byMaterials Science, Institute for Materials

conventional computational fluid dynamics. In

Research (IMR), Tohoku University for the use

particular, we were able to analyze cavitation inof their supercomputers (Project No.20S0026

complex fluids, which was one of the mainfor MASAMUNE-IMR in IMR).

JP19H05718).

We

acknowledgetheobjectives at the beginning of our research. This
achievement would not have been possible

References

without using the ISSP supercomputer system B[1] Y. Asano, H. Watanabe, and H. Noguchi, J.

(Ohtaka).

Chem. Phys. 148 (2018) 144901.

Compared to conventional fluid analysis[2] Y. Asano, H. Watanabe, and H. Noguchi, J.

methods, MD simulation does not require

Chem. Phys. 152 (2020) 034501.

phenomenological models such as equations of[3] Y. Asano, H. Watanabe, and H. Noguchi, J.

state. The advantage of MD simulations is that

Chem. Phys. 153 (2020) 124504.

all physical quantities and phenomena are[4] Y. Asano, H. Watanabe, and H. Noguchi, J.

obtained from molecular motion. Therefore, theChem. Phys. 155 (2021) 014905.

MD simulation is a suitable tool to analyze[5] G. E. Gadd, Nature 206 (1965) 463.

complex flow phenomena, such as the multi-[6] R. E. Arndt, Ann. Rev. Fluid Mech. 34 (2002)

phase and multi-component flow around a143.

propeller.

[7] S. Succi, The Lattice Boltzmann Equation:
For Fluid Dynamics and Beyond, Oxford

Acknowledgment

University Press Inc., New York (2001).

We thank Toshihiro Kawakatsu, Youhei[8] R. Kapral, Adv. Chem. Phys. 140 (2008) 89.

45

Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[9] G. Gompper, T. Ihle, D. M. Kroll, and R. G.
Winkler, Adv. Polym. Sci. 221 (2009) 1.
[10] T. Murashima, ISSP Supercomputer
Activity Report 2015 (2015) 35.
[11] X. Xu and P. Yu, J. Non-Newtonian Fluid
Mech. 229 (2016) 27.
[12] J. R. Cressman, Q. Bailey, and W. I.
Goldburg,

Phys.

Fluids13(2001)

867.

[13] Y. Asano, H. Watanabe, and H. Noguchi, J.
Phys. Soc. Jpn. 88, 075003 (2019).
[14] C. E. Brennen, Cavitation and Bubble
Dynamics, Oxford University Press, (1995).

46