8

Numerical simulation of air/water multiphase flows
for ceramic sanitary ware design by multiple GPUs
A Hybrid Quantum-Classical Simulation
study on the Li Diffusion
in Li-Graphite Intercalation Compounds
A Large-scale Simulation on CFD
in Construction Industry
Estimation of strong ground motion anddamage of indoor objects during giant earthquake

Numerical simulation of air/water
multiphase flows for ceramic sanitaryware design by multiple GPUs
Akio Ikebata* Shinya Yoshida* Feng Xiao***TOTO LTD., Production technology center** Tokyo Institute of Technology, Department of energy sciences

We have been developing an in-house CAE air/water two-phase numerical code for various purposes in design andmanufacturing of plumbing products such as ceramic sanitary wares. In order to re-produce the complex interfacial
flows of air and water with adequate accuracy, large scale computations are required with reliable numerical model,
which is of great challenge. To this end, we have made efforts to improve the numerical schemes and port the code
to the GPU platforms to accelerate the large scale computations for real-case applications. We have implemented
large-scale simulation on the TSUBAME2.0 supercomputer by making effective use of the GPGPU architecture, andachieved significant improvement in both computational performance and simulation results.

Introduction

1

order accuracy and can be combined with the semi-Lagrangian
scheme to get computational efficiency. However the original
CIP scheme doesn’t ensure numerical conser vation, which
may turn out to be a serious problem in many applications,such as interfacial multi-phase flows where mass (or volume)

Being a world-wide leading company, TOTO manufac turesconservation is a basic requirement to resolve small droplets,and provides a wide spectrum of household products, such as

bubbles and water films.

ceramic sanitary wares, bathroom and kitchen fittings. Facing

For the purpose of numerical simulations for real-casethe increasing demands for saving energy and water, TOTO has

product design, we have devised a conservative variant of thealways targeted the high standard in this respect for the design

CIP method, so-called UTI-VSIAM3 scheme [4], and developed an

and manufacture processes. For this purpose, development andaccurate and efficient numerical code for two-phase interfacial

effective use of accurate and reliable CFD (computational fluidm u l t i p h a s e f l o w s . T h e co d e h a s b e e n p o r te d to G P G PU

dynamics) technology as a CAE design tool become extremely

hardware, and tuned on the TSUBAME 2.0 supercomputer fornecessary.

large scale high-resolution simulations.

In the numerical models which are based on Navier-This paper presents our recent ef for ts to improve

Stokes equations and explicitly resolve the moving interfaces

the numerical schemes for multi-phase interfacial flows and toof gas/liquid multiphase flows, like those involving gas bubbles

implement the code on GPGPU architecture.

and liquid films, the advection scheme plays a crucial role forgetting stable and accurate numerical solutions.
MUSCL scheme [10] is perhaps the most representative
a d ve c tio n s ch em e, an d has b e en w id el y a d o pte d in th e

The numerical scheme for advectioncommercial CFD codes. MUSCL scheme is essentially based2on the linear reconstruction and makes use of upwinding andlimiters to get stable and monotonic numerical solutions. As we

We consider the advection equation

will show later, MUSCL scheme appears to be too diffusive toresolve the fine structures which always associate large jumps
and discontinuities in flow fields. Moreover, our experiences
sh ow t h at MUS CL s c h e m e usu a l l y f a i l e d i n re p ro d u c i n greasonable results in violent liquid splashing problems where

where φ denotes the advected quantity, and u,v,w, the components

the MUSCL scheme seems to be more grid dependent even theof the velocity field.
As presented in [5] , a one dimensional conservative CIP

physical problem is geometrically isotropic.
As a completely different alternative, the CIP

02

[1-3]

schemescheme, CIP-CSL (conservative semi-Lagrangian) scheme, uses

is a highly accurate advection scheme which is par ticularly

at least more than two kinds of moments as the computational

more robust and less grid-dependent. CIP scheme uses multi-

variables. In order to make the numerical model computationallymoment concept and updates not only the point value but

efficient, we construct the multidimensional formulation based on

also the gradients of the physical fields. CIP scheme has third-

two kinds of moments [6], i.e. VIA (Volume Integrated Average) andSIA (Surface Integrated Average) as shown in Fig.1 for 3-dimensional

Step 2：model. The one dimensional conservative CIP schemes can beused as the building block through a splitting procedure. So, areconstruction of either quadratic polynomial or rational function
can be build over a single cell. This local reconstruction withadequate accuracy is very important to resolve complex flow

The one-dimensional scheme is firstly used to solve

structures in presence of discontinuities, such as multi-phase flows.

equation (2) in x direction to get the intermediate value for

In practice, the one-dimensional scheme is consistently applied

SIA-X, the SIA for the surface element normal to x direction. The

to SIA and VIA moments in each direction respectively. The SIA isintermediate value of SIA-X is then used to compute the spatial

updated by a semi-Lagrangian solution, while the VIA is computed

discretizations for equation (3) where a simple central difference

from a finite volume formulation of flux form that guarantees

is applied. The SIA-X is then updated to the new time step. Thisrigorously numerical conservation. This formulation, so-called

procedure applies to SIA-Y and SIA-Z. Given the SIA values on

VSIAM3 (Volume/Surface Integrated Average based Multi-Moment

all surface elements, the VIA can be updated by a finite volume

Method), has been further extended to be a general solver for fluid

formulation. The scheme is called UTI-VSIAM3 (Unsplit Time[6]

dynamics .

Integration VSIAM3), which essentially ensures that the updatingof the SIAs in different directions are conducted in a symmetricalmanner [4].

Figure 1 The definition of VIA and SIA moments

In the original VSIAM3 formulation, the one-dimensionalbuilding block is used in a straightforward way for multi-dimensional

Figure 2 Multi-dimensional updating of the SIA moment.

calculation. Alternatively, a multi-dimensional updating canconstructed as follows.

We evaluated the advection schemes by transporting a “T” shape

Shown in Fig.2, we consider the advection computation

initial distribution via a pure anti-clock rotational velocity field.

of SIA-X. The advection equation is divided equivalently into next

The numerical results after one revolution are displayed in Fig.3.

two-step form.

The results from both original VSIAM3 and UTI-VSIAM3 are much

Step 1：superior to the MUSCL scheme. For this pure advection test, there
is not visible difference in the numerical solutions between theoriginal VSIAM3 and UTI-VSIAM3 though, the multi-dimensional
updating in UTI-VSIAM3 can be expected beneficial in real-caseapplications.

03

Numerical simulation of air/water multiphase flows
for ceramic sanitary ware design by multiple GPUs

STAA (Surface Tracking by Artificial Anti-diffusion) [7] scheme.
STAA scheme combines the conservative CIP (CIP-CSL) schemewith an anti-diffusion post-processing step to the VOF function.
The CIP-CSL scheme is able to compute the advection transportwith adequate fidelit y. However, same as other high orderschemes based on finite grid resolutions, the advected field is
usually smeared out due to numerical diffusion. To maintain thesteepness of the VOF function through the whole computation
procedure, correc tion is added to the VOF field af ter theadvection step. To this end, a level set function is generated
from the VOF field, which then gives the normal direction of theinterface. The numerical diffusion from the advection calculationis corrected through an artificial compression in the normal
direction of the inter face, which improves the geometricalfaithfulness of the moving inter face. The STA A scheme isefficient and very easy to use. Fig.4 shows the numerical results
for the aforementioned advection test. The STAA scheme gives
much superior result than the PLIC VOF scheme, especially in
the geometrical shape of the transported field. The thicknessof transitional layer in the STAA calculation is adjustable forpractical use. Our experience in real-case applications shows
Figure 3 Advection benchmark tests.

Interface tracing methodthat the STAA scheme is very stable, robust and reliable.

3

Vo lum e f r a c ti o n is a well sui te d v ar iab l e to i d ent i f y t h einterface in multiphase flow simulation. The VOF (Volume of
Fluid) method uses the volume frac tion ( VOF func tion) for
capturing moving inter faces. Dif ferent from other par ticlebased front tracking methods, VOF method is able to directlycompute complex interfaces of arbitrary topological changes,
such as breaking and merging of the interfaces. Moreover, VOF
schemes automatically guarantees the conservation of the VOF
function, which is of great importance for real-case applications
and makes the VOF among the most popular methods forinterface capturing/tracking. The conventional VOF schemes
adopt geometrical reconstructions which are required then in
estimating the numerical fluxes to update the VOF function.
Usually a PLIC type reconstruction [9] is needed to give adequate
numerical accuracy. However, the PLIC algorithm is complicatedand not easy to implement.
In order to calculate moving interfaces more accurately
and efficiently, we have devised an interface capturing scheme,

04

Figure 4 Interface transport tests.

Porting to GPGPU architecture4

Our numerical code is based on a frac tional step solution
pro ce dure. Following the numeric al s teps for adve c tion,diffusion, gravity, surface tension and other physical processes,a pressure projection step is used to update the velocity field
to the new time level. All computations can be carried out inan explicit manner except the Poisson equation for pressure
projection. The whole code is organized so as to be suitablefor vector processing. Thus, there is not substantial barrier toFigure 5 Comparison of the execution time of 		
PCG Poisson solver on CPU and GPU.

run the code on a parallel environment. In order to implement
the code on parallel machines with distributed memories and
GPGPU accelerators, we have re-written the code thoroughlywith MPI and OpenCL libraries. Although our current target GPU
is nVIDIA CUDA architecture, we make use of OpenCL to allow
the code flexibility to be portable across a variety of GPUs andCPUs, like AMD and Intel’s or other new devices in near future.
It is known that OpenCL usually results in lower performanceabout 10% compared to CUDA. Improvement has been achievedby compiling the code with CUDA Driver API conversion toutilize the native CUDA functions. Furthermore, we have definedsome wrapper functions to simplify the OpenCL instructions in
the fluid code.
Fig.5 plots the comparison of the elapse time between
CPU and GPU. Double precision is used to guarantee the iteration
converging with an adequate accuracy. 8 times speed-up is

Figure 6 The simulation results of the multiphaseflows containing gas, liquid and solid.

achieved on single core and nearly 10 times speed-up on 4 cores.
Fig.6 shows an example of simulating fluid flow in aceramic sanitary ware. The water flowing on the solid surface ofcomplex geometry, as well as the interaction with floating bodies,

Implementation on TSUBAME2.0

was reasonably reproduced. The floating bodies are represented by5

the volume fraction of solid, similar to that used for distinguishingair/water surface[8]. Again, the level set function is used to compute

Since water-saving increasingly becomes a key requirement

the distances among the solid bodies, which is then used in thefor all plumbing fixtures, how to design products that use lesscalculations of contact and rebounce model. All of these are codedwater and have high washing efficiency is the central task we arein a completely portable way to any parallel GPGPU architectures.

facing currently. From the CAE point, we need to resolve the thin

So, we can carry out directly fluid-solid interaction simulations

water film and small bubbles in numerical simulations, which(including gas, liquid and solid) of a relatively large grid size evencan be only realized with finer grid resolution. For example, a

on a personal computer with low-cost GPU cards, which wasmesh-spacing less than 1mm is required to resolve a film withthought impossible unless a supercomputer would be available

2mm thickness. Thus, over ten million mesh elements are usually

before. Establishing such a system is particularly important for theneeded for the whole system, which is quite challenging forpractical use in product-design process, where numbers of casea PC-based hardware. So, we have implemented our code on

studies can be conducted simultaneously on several PCs within

TSBAME2.0 supercomputer in Tokyo Institute of Technology for

acceptable computational time.

large scale simulations of multiphase interface flows.

05

Numerical simulation of air/water multiphase flows
for ceramic sanitary ware design by multiple GPUs

As the computation performance of a single GPU isFig.9 shows the parallel speedup of the Poisson solver

significantly superior to a single CPU, the data communication

for real-case application where 150 million mesh elements

overhead in parallel processing becomes the major barrier towere used. 80 times speedup is obtained on 100 GPUs and 100

get high parallel efficiency. Because the flow path structure

times on 150 GPUs, which shows a reasonable parallel efficiencyof the ceramic ware is ver y complex, it is not easy to getof the present code. Since the Poisson solver is the mosta satisfac tor y parallel scalabilit y if the parallel par tition iscomputationally expensive part in the numerical model for fluid

decomposed along the flow path. On the other hand, a regular

dynamics, the parallel performance shown above reveals that

3D Car tesian decomposition will include a large number ofthe parallelization of the code for GPU architectures are quite

elements that are not actually used in the flow computation. We

promising for large scale simulations.

have finally adopted a flexible 3D partition as shown in Fig.7. Allflow path elements are decomposed along z, y and x directions
automatically, which allows the computational load to be evenly
assigned to each process. Consequently, the halo regions might
be not exactly matched among the neighboring domains, which
complicates the data transfer for parallel processing. We use
asynchronous communication for data transfer based on a predetermined send/receive route. To save the wall-clock time, we
overlap data communication with computation as shown in
Fig.8. All data transfers are finally synchronized by MPI_Waitall
function.

Figure 7 Three dimensional parallel partition.
(conceptual rendering)

Figure 8 Overlapping data communication.

06

Fig.10 displays a real-case simulation of water flows in asanitary ware. All results were also visualized on TSUBAME2.0 with
POV-Ray. Being in agreement with the real design conditions, thewater initially stored in the tank was released with the amountcontrolled by the drain valve. The water flushing through thewhole flow path was simulated. The fine liquid structures andthe entrainment of massive bubbles were adequately captured.
The large scale simulation on TSUBAME2.0 allows us to explicitly
resolve the bubble, thin film and spray of millimeter size, which
drastically improves the accuracy and reliability of the numericalsimulations for the flushing process of sanitary set.

Figure 9 Speedup of PCG solver by multi-GPUs.

Figure 10 Real case simulation of sanitary set
by TSUBAME 2.0 supercomputer.

07

Numerical simulation of air/water multiphase flows
for ceramic sanitary ware design by multiple GPUs

Summary

6[8]

F. Xiao: A Computational Model for Suspended Large Rigid
Body in 3D Unsteady Viscous Flow, J. Comp. Phys., Vol. 155,
pp.348 (1999)[9]

W.J. Rider, D.B. Kothe: Reconstructing Volume Tracking, J.
Comp. Phys., Vol. 141, pp.112 (1998)

For accurate and efficient numerical simulations of air/water[10] B.V. Leer: Towards the Ultimate Conser vative Dif ference

multiphase flows for sanitary ware design, we have developed

Scheme Part V : A Second Order Sequel to Godunov's Method, J.

a reliable numerical model based on the UTI-VSIAM3 and STAA

Comp. Phys., Vol. 32, pp.101 (1979)

schemes. The code has been parallelized by optimizing the loadbalance for the effective flow path region via a flexible threedimensional domain decomposition, and ported to hardware
with GPU accelerators. The implementation of the present
code on the TSUBAME2.0 supercomputer succeeded in directly
resolving the fine flow structures, like bubble, spray and film of
millimeter size in real-case applications. Further efforts will be
devoted to improve the efficiency of the Poisson solver and tophysical models, such as the wetting effects of the solid surface.
More real-case applications will be conduc ted for productdesign.
References[1]

T. Yabe, T. Aoki: A Universal Solver for Hyperbolic-Equations
by Cubic-Polynomial Interpolation 1. One-Dimensional
Solver, Comp. Phys. Comm., Vol. 66, pp.219 (1991)[2]

T. Yab e, T. Ishik aw a, P. Y. Wan g, T. Ao k i, Y. K a d ot a, F.
Ikeda: A Universal Solver for Hyperbolic-Equations by
Cubic- Polynomial Interpolation 2. 2- Dimensional and3-Dimensional Solvers，Comp. Phys. Comm., Vol. 66, pp233(1991)[3]

T. Ao k i: M u l t i - d i m e nsi o n a l Ad ve c t i o n o f CIP (Cu b i cInterpolated Propagation) Scheme, CFD Journal, Vol. 4,
pp.279 (1995)[4]

A. Ikebata, Y. Muraoka, F. Xiao: Multiphase Fluid Simulations
on a Multiple GPGPU PC Using Unsplit Time Integration
VSIAM3, Progress in Nuclear Science and Technology, Vol.2,
pp.491 (2011)[5] T. Yabe, R. Tanaka, T Nakamura, F. Xiao: Exactly Conservative
Semi-Lagrangian Scheme (CIP-CSL) in One Dimension, Mon.
Wea. Rev., Vol. 129, pp.332 (2001)[6]

F. Xiao, A. Ikebata, T. Hasegawa: Numerical Simulations of
Free-interface Fluids by a Multi-integrated Moment Method,
Computers and Structures, Vol. 83, pp.409 (2005)[7] A. Ikebata, F. Xiao: Development of Conser vative Front
Capturing Scheme and Applications to Multi-f luid
Simulations, Proc. Annual Conference of Japan Society of
Mechanical Engineering, 3, pp.301 (2002) (in Japanese).

08

A Hybrid Quantum-Classical Simulation
study on the Li Diffusion in Li-Graphite
Intercalation Compounds
Nobuko Ohba* Shuji Ogata** Takahisa Kouno*** Ryoji Asahi** Toyota Central R&D Labs., Inc.

** Nagoya Institute of Technology*** Institute for Solid State Physics; The University of Tokyo

The Li diffusion process in the Li-graphite intercalation compound (Li-GIC) is investigated by the hybrid quantum(QM)-classical (CL) simulation. The region that includes the inserted Li and neighboring C atoms is treated withelectronic state by the QM calculation and then embedded into the CL system of graphite described with anempirical interatomic potential model. A series of the hybrid QM-CL simulation runs on the dynamics of a Li-ion in
the Li-GIC at constant temperature for various values of the averaged inter-layer distance of graphite is performed.
We thereby find that the Li diffusivity is suppressed substantially when the inter-layer distance is compressed bya few percent from the equilibrium value. On the other hand, in the equilibrium and stretched cases, the diffusive
motion of the Li-ion is composed of ballistic and hopping modes. It is found that the Li-ion existing around the
middle of the upper and lower C-layers diffuses fast. Therefore as a result of trying to control the Li position between
neighboring two layers of graphite using an external electric field, we find that in-plane diffusivity of the Li-ion isenhanced by suitable amplitude and frequency of the electric field perpendicular to the C-layers.

Introduction

1(i) The hybrid method can investigate the diffusion process of
the Li-GIC with the influence of surrounding C atoms taken
into consideration at a reasonable computational cost. On
the other hand, the CL molecular dynamics (MD) method
cannot describe generally the chemical reactions such as

The graphite can form various intercalation compounds by

the charge transfer between the Li and C atoms. Also theincorporating atoms (ions) and small molecules between thefirst-principles MD method that uses the DFT for electronic

carbon (C) layers. Among these, the Li-graphite intercalation

structure calculation cannot treat a large-scale graphite

co mp o un d (L i - GI C ) is p u t to p r a c ti c al us e as a n e gati vesystem required to study the migration of the Li-ion.

electrode of the Li-ion rechargeable battery. Among the key(ii) The dispersion forces for the inter-layer interac tion ofprocesses, the transport process of the Li-ions in the graphite

the graphite can be incorporated easily through the CL

layers is related directly to the power performance of the Li-interatomic potential model in the hybrid QM-CL method.

ion batteries. With the Li insertion process, the Li-ion creates

Remember that the dispersion force cannot be taken intoa long-ranged stress field around itself by expanding theconsideration in the conventional DFT method.

inter-layer distance of the graphite by approximately 10 % . To(iii) It is not necessary to construct an empirical interatomic

take into account such a long-ranged stress field in the first-

potential between the Li and C atoms. If we succeed toprinciples simulation of the Li diffusion, the hybrid quantum

construct the proper interatomic potential between the(QM)-classical (CL) simulation code is developed.

Li and C, we can simulate the Li diffusion in the Li-GIC by

For understanding the stress dependence of the Li

the CL-MD method. However, as we will show in § 3.2, thediffusivity in the Li-GIC, we firstly perform a series of the hybrid

activation energy for the Li diffusion is quite small. It is not

QM-CL simulation runs for the dynamics of a single Li-ion ineasy to construct a precise potential between the Li and C

the graphite at temperature of 423 K for various values of theatoms for our purposes.

averaged inter-layer distance. In the hybrid QM-CL scheme,

More over, we will analy ze the relation b et we en

the reaction region where the electronic structure should be

the vertical position and diffusivity of the Li-ion in graphite.

treated by a highly accurate calculation technique such as the

Motivate d by the result s of the present analyses, we will

density functional theory (DFT), is embedded in a CL system

perform additional hybrid simulations to demonstrate enhancedof atoms described with an empirical interaction model. It isthermal diffusivity of the Li-ion by applying various settings

expected that we can simulate the realistic large-scale systemof the amplitude and frequency of alternating electric fields

to reproduce the physical phenomenon of our interest by using

perpendicular to C-layers.

the hybrid method. Remarkable merits of using the hybrid QMCL method for the present simulation of the Li-ion in the Li-GIC
are the following.

09

A Hybrid Quantum-Classical Simulation study
on the Li Diffusion in Li-Graphite Intercalation Compounds

2

Methodnodes. The multi-grid method [9] is employed for accelerationof the convergence of the long-wavelength components ofthe data on the mesh points. The grid spacing h=0.45 a.u. (1 a.u.
≈ 0.5292 Å ), which corresponds to the cutoff energy ( π/h ) 2 ～49 Ry(1 Ry = 13.6 eV ) in the plane waves DFT method. In addition,

2.1 Buffered cluster method

the smaller grid spacing of h/3 is used only around the atoms.

In the hybrid QM-CL simulation method, a cluster of atoms

external electric field E, we consider the following Kohn-Shamis selected from the total system as the QM region. Artificial

equation in the DFT:

For the QM re gion with buf fer atoms under the

dangling bonds therefore form at the QM-CL boundary. Possible
influence of the dangling bonds on the electronic states and theatomic forces should be minimized. The link-atom method thatuses the hydrogen atoms for termination of the dangling bonds

with an eigen orbital and energy of

is often adopted to couple the QM and CL regions. Depending

atomic unit is used). Here, v eff is the effective potential energy

on the selection of the QM region, however, using the link atom

for an electron at r without the external electric field. Charge

method results in significant relaxation of the artificial surface ofneutrality is assumed in the calculation. The force on atom- I is

the QM region and unwanted deformation of the whole system

obtained by adding F I = Z I E ( Z I is the ion charge) to the atomic

from the original configuration. We therefore adopt the buffered

force caused by electron-electron, electron-ion, and ion-ion

cluster method (BCM) [1] , which requires no link-atom and isinteractions.

insensitive to the selection of the QM region. To use the BCM, we
put additional atoms of either H or C (i.e., constituent species),, respectively (the2.3 CL calculation

called the buffer atoms, at the QM-CL boundaries to terminate theIn the hybrid QM-CL simulation method, the CL-MD calculation

dangling bonds. The positions of the buffer atoms are determinedis performed for the CL region, that is, the whole system minus

so as to minimize the potential energy of the corresponding

the QM region. The velocity Verlet algorithm is used to integrate

cluster in CL calculation under the constraint of fixed the positions

the New ton’s e quations of motion for all the atoms. T heof atoms in the QM region. Note that the positions of the buffer

Brenner’s interatomic potential [10] is adopted for the C atoms ofatoms are not relaxed in the QM calculation. Thereby artificial

the graphite. We employ the parameter values listed in Table I

surface relaxation of the QM region is suppressed in the BCM.

in Ref. [10] except for the equilibrium distance. The equilibrium
distance is set to 1.33116 Å so that the calculated equilibrium

2.2 QM calculation

lat tice constant of the graphite becomes equal to the one

For the QM calculation in the hybrid QM-CL simulation, we use

calculated with the QM calculation, which is slightly longer than

the real-space grid based DFT (RGDFT) method in which thethat in Ref. [10]. Such a fine tuning of the parameters in the CL

Laplacian operations in the Kohn-Sham and Poisson equations

potential is necessary in order not to create artificial stress atare evaluated with the finite difference method. The RGDFT is

the QM-CL boundary in the hybrid QM-CL simulation method.

well-suited to the present setting of the free boundary condition

In the Brenner-type potential, the cut-off distance ofin the QM calculation and is highly parallelizable. The details of

the interaction is 2.0 Å. It means that only the interaction between

the algorithm were described in Refs. [2-4]. The Troullier-Martin-

the nearest neighbor atoms belonging to the same layer istype normconserving pseudopotentials [5] are used to describe

considered. To take into account the inter-layer interaction, we

the interaction of the valence electrons and nuclei (ions). The

construct the interatomic potential model V vdw [11] based thegeneralized gradient approximation formula with partial core

Lennard-Jones potential for the dispersion (or the van der Waals)

correction [6] introduced by Perdew, Burke, and Ernzerhof [7] isforces for C atoms in different layers to add it to the CL potential.

adopted to the exchange-correlation energy term. The KohnSham orbitals and Hartree field are represented on the uniform
Car tesian mesh points. The four th- order finite - dif ference[8]

method is used for the Laplacian operation .

10

and

2.4 Calculation model
The simulation system is composed of 3072 C atoms and a single
Li-ion under the periodic boundary condition as shown in Fig.1.

For parallel computation, those data on the mesh

The graphite assumes the AB sequence since we are interested

points are divided into domains to be stored in the computein the situation of low Li-density. The x - and y -axes are set alongthe layer, while the z -axis perpendicular to the layer. The hybrid
QM-CL simulation at T=0 K gives the simulation system of L x

Stress dependence of the diffusivity of=34.44 Å, L y = 29.83 Å, and L z = 26.97 Å as the equilibrium values.

the Li in graphite

Since the thermal expansion along the layer is quite small,3

we use the same values of (L x , L y ) at finite temperatures also.
On the other hand, the thermal expansion along the z-axis is

3.1 Li diffusion in changing interlayer space of graphite

substantial. It is provided that the coefficient of linear thermal
expansion of the pure graphite in z-direction is 6.7 × 10 K − 1 by

The trajectory of the Li-ion during 10 ps in each run is plotted in

the separate CL-MD simulation with the present CL potential.

Fig.2; the initial and final positions of the Li ion are depicted withWe assume that the same expansion coefficient can be applied

the open arrow and the large red sphere, respectively. The green

to the present system since the Li density is too small to affectspheres are the C atoms in the QM region. When the graphite

the system size. Therefore, multiplying the linear expansion

layers assume the equilibrium distance or stretched one, the Li

factor for 423 K to L z = 26.97 Å for T = 0, we find the equilibrium

ion migrates over a wide area. In the equilibrium and stretchedvalue of L z = 27.74 Å for T = 423 K. We investigate the Li diffusioncases (0.6 % and 2.9 % ) in Fig.2, we observe both hopping andprocess when changes in L z the original value of 27.74 Å are -3.9 ,ballistic modes of the Li ion motion. In the hopping mode, the Li-1.7 % (compressed) , 0.6 % (equilibrium) and 2.9 % (stretched).

ion moves to avoid the places at which two C atoms belonging-5

to different layers assume the same x-y positions (see, Fig.4).
However, as remarked in Fig.2 with black arrows, we find the ion
in the ballistic mode can pass through such a place. Detailed
analyses of the ballistic-motion events show that the inter-layer
distances at near places of the Li-ion toward the direction of its
motion are always larger than 3.85 Å. We, in fact, find that thethermal fluctuation produces such a local stretching frequently
at various places in those simulation cases. On the other hand,in the compressed case the Li-ion diffuses in the hopping mode
only and appears to be confined in the deformed layers (i.e., thecage effect).
Figure 3 shows the mean square displacements (MSD)of the Li-ion as functions of time in the four cases (-3.9 % ,-1.7 % ,0.6 % , and 2.9 % changes of L z ) at T = 423 K. Since there exists only a

Figure 1 The side (x−y ) view of the total simulation
system. The large red and medium green
spheres are respectively the Li and C atoms
in the QM region. The small blue spheres
are the C atoms in the CL region.

single Li-ion in the system, we shift the time origin of the MSD
by ever y 0.01ps to increase the number of data samples forbetter statistics. The total simulation time is 12ps for Li each run;the time step is 1.0 fs. The curves in Fig.3 are not linear enoughto estimate the diffusion coefficients at long times. However,characteristic dependence of the Li diffusivity on the changeof L z is observed as follows. The MSD’s in the equilibrium (0.6 % )and stretched ( 2.9 % ) cases are larger than that in the compressed(-1.7 % and -3.9 % ) cases. However, the long-time diffusion in 2.9 %of stretched case tends to become slightly slow. This is because
the Li-ion is trapped near the upper or lower graphite layer
and the details will be provided later. In a compressed case, theclear change of the MSD from the linear to saturation behavior
is obser ved. The mechanism of the saturation behavior is
explained in terms of the cage effect.

11

A Hybrid Quantum-Classical Simulation study
on the Li Diffusion in Li-Graphite Intercalation Compounds

Figure 4 The Li ion at the ground and transition
states in the hopping diffusion in thegraphite. The Li ion is sandwiched by
two C-layers.
Figure 2 The trajectory of the Li ion during 10 psviewed from z-direction, obtained in the
hybrid QM-CL simulation. The four casesof the change of L z are considered. The
initial and final position of the Li ion isdepicted by the open arrow and the large
red sphere, respectively. The green spheres
are the C atoms in the QM region. The blackarrows depict the places where the Li ion
passes through the places at which two C
atoms belonging to different layers assume
the same x- y positions.

3.2 Activation energy
In order to understand the mechanisms of the Li dif fusivity,
the activation energy for the hopping diffusion of the Li ion isevaluated using the hybrid QM-CL simulation method. Figure
4 depicts the positions of the Li ion in the ground state (GS)and the transition state (TS). The activation energy, which is
calculated by subtracting the energy in the GS from that in TS,
in each case of the L z change is plotted in Fig.5. The activationenergy obtained from the total energy, i.e., the filled sphere in
Fig.5, includes both contributions of the QM and CL regions. We
find that the activation energy is quite small, less than 0.07eV,
irrespective of L z . Therefore, the Li ion can diffuse easily by thethermal fluctuation. As the inter-layer distance becomes shorter,
the activation energy estimated from the total energy increases
subs t antiall y, w hich sup p o r t s o ur f in din g o f sup p ress e ddiffusivity by compression. The contribution of the QM region(i.e., around the Li ion) to the activation energy changes little
when the system is compressed or stretched. On the other hand,the contribution of the CL region (far from the Li ion) increases
as the system is compressed, which means that the distortion of
the C atoms far from the Li ion affects the activation energy.
To analyze possible error in the activation energy, we
calculate the energy with a different method in a similar setting.

Figure 3 The mean square displacements of the
Li-ion as functions of time in the hybrid
QM-CL simulation at 423 K for thecompressed (-3.9 % and -1.7 %), 		
equilibrium (0.6 %) and stretched (2.9 %)values of L z .

We use the projector augmented wave method for the DFT
calculation implemented in the VASP code [12], which is expected
to have higher accuracy since the core electrons are treated also
in addition to the valence electrons. We prepare the supercell, 3
×3× 2 of the unit cell, of the graphite in the AB-sequence. Thena single Li atom is inserted. Hence the total system contains 72 C

12

atoms and 1 Li atom. The calculated results by the VASP are listed
in Table 1 for both GS and TS. The activation energy evaluatedEnhancement of Li diffusion by

from the electronic energy is 0.09 eV, which is the same order

external electric fieldobtained by the hybrid QM-CL simulation method at zero strain.

4

Considering the difference in the target system used in the two
methods, we state that no substantial difference in the activation

4.1 The relation between the vertical position and diffusionenergy exists between the two results. In addition we calculate

behavior of Li at E=0

the activation free-energy using the phonon package MedeA [13].
Relatively large contribution of the phonon entropy TΔS = - 0.14 eVFirstly, we analyze the relation between the vertical position andis found at 400 K. Finally the difference in the Gibbs free-energy

diffusion behavior of the Li ion observed in the hybrid simulation

between the GS and TS becomes 0.04 eV.

run with a zero electric field, i.e., E=0 . The stretched case (2.9 % ,L z =28.55 Å) in which the Li-ion diffuses well is adopted. Figure 6(middle) shows the time evolutions of the distances d u and d lbetween the Li ion and the center-of-mass positions of QM-region
C atoms in the upper and lower layers, respectively. The time
evolutions of absolute values of Li displacement in the x-, y-, andz -directions |(R Li(t) − R Li(0))x,y,z | are shown in Fig.6 (bottom). In most of

the period of t=0 to 1 ps, d l is shorter than d u, which means that
the Li-ion is trapped near the lower layer. During the period, theLi-ion hardly diffuses, as shown in Fig.6 (bottom). In the following
period of t=1.0 to 1.2 ps, on the other hand, d l and d u are nearly
equal and the Li dif fusivity is relatedly high. Such a relation
between d l , d u and the in-plane displacements of the Li-ion can
Figure 5 The activation energy for the hopping diffusionof the Li ion in the graphite at T=0K calculatedusing the hybrid QM-CL simulation method.
The filled circles with the solid line represent
the data evaluated using the total energy. The
open circles with the dashed line represent thedata using the energy of the QM region only.
The strain corresponds to the variation of L z
from the equilibrium distance of 26.97Åat 0K.
The error bars indicate the estimated overall
errors relating to the convergence fluctuationin the real-space DFT method.

Table 1

The results of various quantities at the groundstate (GS) and transition state (TS) calculated
400K
400K
using the VASP code. The E ele, E vib, S vib , and
400K
400K
400K
G
= E ele +E vib +pV −TS vib are the electronicenergy at T=0K, the vibrational energy including
the zero-point energy, the phonon entropy, and
the Gibbs free-energy at T=400 K, respectively.
The electronic energy is measured in reference
to the standard states of the atomic elements.

also be seen at approximately t=1.6 ps and 4.2 ps in Fig.6. We maystate that Li diffusivity is enhanced when the Li-ion escapes fromone of the two stable vertical sites to move around the middle ofthe upper and lower layers where a resisting force on the Li ion forthe perpendicular ( x or y ) motion is relatively weak.

Figure 6 Diffusion behavior of Li-ion in graphite at T=423
K with zero external electric field, observed in
hybrid QM-CL simulation run. The time
evolutions of the external electric field, du anddl, and the absolute values of the Li
displacements in three directions are depicted
in top, middle, and bottom panels, respectively.

13

A Hybrid Quantum-Classical Simulation study
on the Li Diffusion in Li-Graphite Intercalation Compounds

4.2 Li diffusion in applying external electric field

The valence-electron density distribution in the QM
region applied the electric field which of amplitude is 0.43 V/Å

Motivated by the finding stated above, we think of increasing

statically is shown in Fig.9. If the amplitude of electric field is too

the mean probability of the Li residence around the middle

large, valence-electron cloud have possibilities to greatly slip offof the upper and lower C-layers to enhance the diffusivity offrom the atomic nucleus position. Present electric field amplitude

the Li-ion in graphite by applying an alternating electric fieldof 0.43 V/Å is a considerably strong, but the electronic distribution

perpendicular to C-layers.

remains around the carbon atom of the QM region.

Figure 7 shows the MSD of the Li-ion as functions of2

No experimental report exists about the vibrationaltime〈|(R Li(t) − R Li(0))| 〉obtained in the hybrid QM-CL simulation

frequency peak of Li-GIC in the range of 1 ∼ 2 THz. For graphite,runs at T = 423 K with electric field amplitudes E ampz =0.0 , 0.14 ,the low-frequency mode of E 2g(1) is observed at 1.26 THz in theand 0.43 V/Å and frequencies f =0.4, 0.8 , and 1.7 THz. For better

Raman measurement [14] . If the oscillation frequenc y of thestatistics, the time origins are set every 10 fs in a simulation

external electric field is tuned to the mode, it is likely that therun with a total of 5 ps. As E ampz increases with f =0.4 THz, the Livertical vibration of the layers in graphite enhances Li diffusivity

diffusivity gets lower. This is because the Li-ion experiences afurther by increasing the probability of the Li residence aroundforce to stay around one of the two stable vertical sites exerted

the middle of the upper and lower layers. We note that no suchby the vertical electric fields at such a relatively low frequency.

vibration of the layers is observed in the present simulation

From the Fig.5, the activation energy for the hopping diffusion

runs. Technologies for generating high-intensity radiation inof the Li ion is estimated as about 0.04 eV for the present system

the teraher tz range have been advancing rapidly in recent

size of L z = 28.55 Å. Since the vertical distance between the stable

years, for example, using the gyrotron system [15]. Experiments tovertical site and the transition state site is about 0.2 Å, E amp;thresz =demonstrate our present prediction of enhanced Li diffusivity in

0.04/0.2 = 0.2 V/Å may be the threshold value for such a verticalgraphite by using an alternating vertical electric field should betransition of the Li-ion to occur without the thermal energy. The

possible in the near future.

Li-ion is expected to escape quickly from one of the two stablevertical sites following the time evolution of the electric field if
E ampz > E amp;thresz .
We therefore analyze the f dependence of〈|(R Li(t) −R Li(0)) x,y,z |2〉for E ampz = 0.43 V/Å, which is larger than E amp;thres z . As
shown in Fig. 7, the Li diffusivity at f = 0.8 THz is similar to that at
f = 0.4 THz and that the Li diffusivity at f = 1.7 THz is remarkably

high. The time evolutions of the Li-C distances d u and d l forE ampz = 0.43 V/Å at f = 0.4 and 1.7 THz are shown in Fig.8. At f = 0.4
THz, the vertical position of the Li ion is high (low) when theelectric field directs upward (downward) as seen in Fig.8(a).
Since the Li-ion hardly diffuses when it stays at one of the twostable vertical sites, the effective diffusivity becomes low atsuch a low frequency. At a high frequency of f = 1.7 THz, the Liion changes its vertical position in a chaotic manner with thelarge and alternating vertical forces on it owing to the external
electric field, as shown in Fig. 8(b). It appears that the verticalmotion of the Li-ion cannot catch up with such a fast change inthe direction of the external electric field. In other words, theLi ion experiences a retarding force from the external electricfield before arriving at a stable vertical site. The Li ion therefore
increases its probability of residing around the middle of theupper and lower layers, resulting in a significant enhancementof Li diffusivity.

14

Figure 7 Time evolutions of mean squared
displacements of Li ion obtained in hybrid
QM-CL simulation runs at T = 423 K withvarious external electric fields. For better
statistics, the time origins are set every 10
fs in a total simulation period of 5 ps.

(a) At frequency f=0.4 THz

Summary5

We have applied the hybrid QM-CL simulation code to analyze
the stress dependence of the dif fusivity of a single Li-ion in
the graphite at T=423 K by changing the averaged inter-layer
distance (by plus and minus a few percent). Thereby we havefound that the Li diffusivity is suppressed significantly in thecompressed case, while in the stretched case, the Li-ion motion
has shown both hopping and ballistic modes. On the other(b) At frequency f=1.7 THz

hand, in the compressed case, the Li ion has diffused throughthe hopping mode only and has been found to be confined ina relatively small area at longer times; its mechanism has been
explained in terms of the cage effect. From our findings about
the stress-dependence of the Li diffusivity in the graphite, we
think a possible mechanism for the drastic lowering of the Li
diffusivity at particular densities of Li observed experimentally in
the Li-GIC.
We have found that the in-plane diffusivity of the Liion at T = 423 K is enhanced significantly by the electric field
perpendicular to the C-layers if the amplitude E amp z > 0.2 V/Å
and the frequency f is as high as 1.7 THz. The mechanisms of
the enhanced diffusivity have been explained in terms of the

Figure 8 Same as Fig.6, but with the external
electric field amplitude E z =0.43 V/Å .
The right side graphs of the middle panel
show probability distributions of d u and d l,which are obtained by taking average of 5 ps.

increased probability of the Li residence around the middleof the upper and lower layers, resulting in a weak interaction
between the Li ion and the graphite layer.
In the present simulation, the real-space grid based
DFT (RGDFT) code suited to the parallel computer has been
adopted for the QM calculation in hybrid simulation. Fromthe v iew p oint of comput ational cos t , the RGDF T c annot
handle large-scale system including many atoms, because its
computational time increases in proportion to N3 (N: the numberof atoms including in the system). Therefore we have developed
the novel order-N Divide-and-Conquer (DC) type RGDFT code,
in which we divide the target system into several domains [16].
Figure 10 shows a benchmark test of the DC-RGDFT code on
TSUBAME2.0 for large-scale Al cluster. Five calculations for thesystem including the Al atoms of 48, 384, 1296, 3072, and 6000
have been performed with dividing 1, 8, 27, 64, and 125 domains,respectively. The number of SCF iteration is fixed to 25. The

Figure 9 The valence-electron density on a y-planein the QM region applied the electric field
which amplitude of 0.43 V/Å. The black
sphere is the Li ion, while the gray spheres
the C atoms. The density less than 0.001
a.u.－3 is omitted.

computation power increases linearly with N in order to treat 4
atoms per 1 CPU-core. It is confirmed that the total calculation
time using the DC-RGDFT scale linearly with N, while that of
RGDFT increase proportionally with N 3. The use of DC-RGDFT
for larger scale QM region in hybrid QM-CL simulation is now in

15

A Hybrid Quantum-Classical Simulation study
on the Li Diffusion in Li-Graphite Intercalation Compounds

progress. We are working on the case of existing multiple Li’s

[8]dynamics will be clarified.

J. R. Chelikowsky, N. Troullier, and Y. Saad, Phys. Rev. Lett. 72,
1240 (1994).

in the Li-GIC, and the effects the Li-Li interaction to the Li-ion[9]

A. Brandt, Math. Comput. 31, 333 (1977).

[10] D.W. Brenner, Phys. Rev. B 42, 9458 (1990); 46, 1948(E) (1992).
[11] N. Ohba, S. Ogata, T. Tamura, S. Yamakawa and R. Asahi,
CMES-Comp. Model. Eng. 75, 247 (2011).
[12] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
[13] K. Parlinski, Medea-phonon, materials design, inc., 2003,
based on K. Parlinski, phonon 3.11, 2002.
[14] R. Kostic, M. Miric, T. Radic, M. Radovic, R. Gajic, and Z. V.
Popovic: Acta Phys. Pol. A 116, 718 (2009).
[15] T. Saito, T. Notake, Y. Tatematsu, A. Fujii, S. Ogasawara, L.
Agusu, T. Idehara, S. Kubo, T. Shimozuma, K. Tanaka, M.
Nishiura, K. Kawahata, and V. N. Manuilov: J. Phys. Conf. Ser.
227, 012013 (2010).
[16] N. Ohba, S. Ogata, T. Kouno, T. Tamura and R. Kobayashi,
Comput. Phys. Commun. 183, 1664 (2012).

Figure 10 The total calculation time for Al cluster
system of various numbers of atoms N
by the RGDFT and DC-RGDFT method.

Acknowledgements
Some numerical calculations were carried out on the TSUBAME2.0
supercomputer in the Tokyo Institute of Technology supportedby the MEXT Open Advanced Research Facilities Initiative. This
research was supported in part by KAKENHI, Grant-in-Aid for
Scientific Research (B) 233100074 from The Ministry of Education,
Culture, Sports, Science and Technology (MEXT).
References[1]

S. Ogata, Phys. Rev. B 72, 045348 (2005).

[2]

S. Ogata, E. Lidorikis, F. Shimojo, A. Nakano, P. Vashishta, and
R. K. Kalia, Comput. Phys. Commun. 138, 143 (2001).

[3]

S. Ogata, F. Shimojo, R. K. Kalia, A. Nakano, and P. Vashishta,
Comput. Phys. Commun. 149, 30 (2002).

[4]

F. Shimojo, T. J. Campbell, R. K. Kalia, A. Nakano, S. Ogata, P.
Vashishta, and K. Tsuruta, Future Gen. Comp. Sys. 17, 279 (2000).

[5]N. Troullier and J. L. Martins, Phy. Rev. B 43, 1993 (1991).

[6]

S. G. Louie, S. Froyen, and M. L. Cohen, Phy. Rev. B 26, 1738
(1982).

[7] J. P. Perdew, K. Burke, and M. Ernzerhof , Phy. Rev. Lett. 77,
3865 (1996).

16

A Large-scale Simulation on CFD
in Construction Industry
Pham Van Phuc* Tsuyoshi Nozu* Hirotoshi Kikuchi* Kazuki Hibi*** Center for Advanced Computational Engineering, Shimizu Institute of Technology** Numerical Flow Design CO., LTD.

Predicting the performance of computer in a CFD application is an important issue in the practical fields. In thispaper, the performance of a single compute node, GPU and computing networks by year have been investigated.
Large-scale simulation in massively parallel process is necessar y to maximize the performance of computer atpresent. The usefulness of large-scale simulation in using high resolution grid in order to achieve more accuracy
results in comparison with experimental results, in using small grid size to develop a more elaborate numericalmodel reproducing much more flexible shape of object, or to simulate a wide urban area considering the effectsuch as surrounding complex terrain to buildings are also discussed together.

Introduction

1

Trends in performance of computerin a CFD application2

Improving the performance of supercomputers [1] at a rate ofUnderstanding the computational performance of computer

1000 times in the last 10 years has been realized by using theis impor tant in practical applications. In order to clarif y its

modern devices such as GPU and taking advantage of high-per formance, we investigate the changing by year of the

speed and large bandwidth network with many compute nodes.

per formance on single compute node, on the net work of

However, what has been changed or should to be changed incomputer systems, and the scale of numerical model in a

the daily use of computer, the issue has not been discussed

practical CFD application. Because of small different performancemuch in the field of manufacturing and other industries.

between CFD software using the equivalent analytical methods,

In the construction industry, building a safe, secure andcomfortable living environment under considering the effectsof surrounding environment is most important. Many problems

the open source OpenFOAM [2] is used in this study.
2.1 Performance on single compute node [3]such as estimation of wind load, tsunami load on buildings, andFig.1 shows the increasing of the performance of computer by

evaluation of thermal and wind environment around them need

year on a CFD application of wind environmental analysis for an

to be solved together. While the experimental and empiricalurban area with a 10 million unstructured mesh. In the figure, theevaluations are limited, numerical simulation is found as anvertical axis illustrates the speedup has been normalized by theef fective method to solve these problems. Recently, duringcomputed time of one CPU core or a compute node at the year of

the wide spread of BIM (Building Information Modeling) with

2006. In here, the CPU AMD Opteron in 2006 which was adopted

large-scale three-dimensional data, the numerical model itself

in TSUBAME1.0, the Intel CPU Westmere-EP in 2010 which was

becomes much more flexible. Simulating a large-scale model

adopted in TSUBAME2.0, the Intel CPU Xeon in 2008 and the Intelis more challenging and it comes with a high compute cost.

CPU SandyBridge in 2012 are considered. From 2008 onward,

Maximize the computational capacity of computer is a necessary

although the clock rate of CPU had not changed so much, theand urgent issue nowadays.

speedup of a CPU core is now still increased due to the improved

In this paper, the authors try to discern the computational

width of the CPU's data bus, the latency of the memory, theability of computer system in an application of computational fluid

cache architecture and improved floating-point arithmetic unit.

dynamics (CFD) from a practical standpoint. Some applications

But, it is difficult to improve the performance of a single CPU core

are also discussed in here to mark the importance of large-scale

on a serial application in the future because the speedup by year

simulation in our construction industry.

has increased much more slowly. In the other hand, the speedupof a compute node using full cores in a parallel computation
is growing steadily at about 1.5 to 2 times each 2 years with anincrease in the number of cores every year. Therefore, parallel
computation in a compute node is now an important factor to
maximize the computational capacity of computer. However,17

A Large-scale Simulation on CFD
in Construction Industry

the result also indicates that improving the performance of acompute node in 10 years would be less than 32 times.

2.2 Performance on massively parallel computation [3]Fig.2 shows the speedup of GPU to CPU core for the

Recently, because of high performance of computer network,

solvers of a pressure Poisson equation that takes a major part ofusing many compute nodes in massively parallel computation isthe computational time in CFD [4]. Some solvers using GPU can

expected in order to improve the computational performanceobtain high performance in comparison with using a CPU coreon CFD. Fig.3 shows the speedup by node of the Infiniband

in the same time. Therefore, using new devices such as GPU

SDR network in 2006, and the Infiniband QDR network in 2010,

and developing a new appropriate computational method will

which was adopted in TSUBAME1.0 (2 SDR lines with bandwidth

be an important factor in the improvement of computationalof 20GBps) and TSUBAME2.0 (2 QDR lines with bandwidth of

performance on CFD in the future.

80GBps) respectively. Two models with a 10 million and a 400
million unstructured mesh cells, are used. For the 10million
mesh model, the performance is improved by network in 6 times
in 4years from 2006 to 2010. The speedup also grows up much
more with the increase of the number of compute nodes in thelarge model of 400million mesh. Therefore, large-scale modelin massively parallel computation is necessary to obtain a highperformance on CFD application and to maximize the capacityof computer system. Moreover, in considering the performanceof a single compute node as mentioned above, it could be said
that it is possible to improve the performance of computer at
the rate 1000 times in 10 years on a practical CFD application.

Figure 1 Increasing of CPU performance by year ina CFD application

Figure 3 Speedup in different computer networks
and number of meshes of a CFD model
Figure 2 Speedup of GPU by CPU core in a pressure
Poisson equation of a CFD application18

Possibility of large-scale simulation

3

Previous chapter has discussed the computational performanceof CFD focused on the single compute node, network computer
hardware. It has been shown that a large-scale simulation in
massively parallel computation is possible and necessar y atpresent. In this chapter, some practical CFD applications will be
discussed in order to clarify the importance of the large-scale

Figure 4 Computational domain of a wind tunnel

simulation in the construction industry.
A simulation becomes larger when someone wants
to use more high resolution of grid in order to improve morehigh precision of calculated result, to develop a more elaborate
numerical model in order to reproduce a much more flexible
shape of an objec t, or to simulate a wide -area in order toconsider the ef fects of surrounding environment. In here, a
reproducibility of an experimental wind tunnel, an elaboratemesh generation for the buildings with complex 3D shapes andan analysis of three-dimensional tsunami in wide-area using VOF
method are discussed.
3.1 Reproduction of a boundary layer wind tunnel experiment[5]Fig.4 illustrates a computational domain 30m long of a LES
simulation, which was constructed to model conscientiously
a boundary-layer wind tunnel measurement section including

Figure 5 Power spectrum densities of wind flowspires, fences and block roughness. The simulation is carriedout for three models with 1 million mesh cells (model A), 10

generally showed good agreement with the experimental ones

million mesh cells (model B) and 100 million mesh cells (model

for both models and identified the sensitive local pressure zones

C). Fig.5 shows the mesh sizes, the number of mesh cells and theoccurring on roof and wall surfaces. However, it underestimatedcomparison between the experimental and calculated power

the pressure in narrow local pressure zones. Using higher

spectral densities of turbulence wind obtained inside the windresolution mesh, the result becomes better.

tunnel. The result of the model A is significantly lower in the

Moreover, Fig.7 shows the time-averaged streamlines

high frequency region whose dimensionless frequency nLx/

of th e w in d f low p at ter ns aroun d th es e m o dels . For the

UH is less than 1. Using higher grid resolution such as model B

30million mesh model, the flow at ground level separates fromand C, the calculated results become better. Moreover, it also

the windward corner of the lower storey, descending close toshows that the numerical simulation requires approximately 100

the windward wall, and flows through openings beneath thebillion mesh cells in extrapolation in order to fully reproduce theupper storey. However, using 120million refine mesh model, an

experimental results.

upward new cone-like vortex close to the side wall of the upper

In addition, Fig.6 illustrates the distributions of means to rey has b e en fo un d to g eth er. A s a re sul t , larg e - s c al e

and negative peak external pressure coefficients on the roofs

simulation using higher grid resolution is necessary not only

and walls of a typical building model obtained from the pressurein an improvement of the calculation accuracy but also in theexperiment and numerical simulation of t wo models withelucidation of a new mechanism in the flow field.

30milliion mesh and 120million refine mesh using higher grid
resolution. These coefficients are required in the assessment ofwind load on claddings of buildings. The numerical simulation

19

A Large-scale Simulation on CFD
in Construction Industry

3.2 An elaborate mesh generation for thebuildings with complex 3D shapes [3]Using advanced three-dimensional data in the constructionindustry is accelerated in recent years. The complex and detail
3D shapes of buildings are used increasingly in design due tomany different needs. Until now, in order to obtain a mesh of acomputational model of just only several million mesh, which islimited on the computer resource and its computational ability, it
is necessary to use the mesh size of about several ten centimeter
to meter. To generate the mesh, it is usual to clean and simplify
the complex 3D shape data that suitable for the CFD. This work
is often difficult and takes a lot of time depending on the detailsof building shapes.
It should be known that mesh generation becomeseasier when mesh size is much smaller than the scale of thetarget model. Fig.8 shows an example of complex 3D buildings.
Several centimeter of mesh size has been adopted in the meshgeneration process. The mesh is generated automatically
including both the indoor and outdoor sides as shown in Fg.9.
Although the number of mesh reaches to 200milion cells, to bea large-scale model, building model in detail such as handrails
Figure 6 Pressure distributions on a building

and interior are reproduced as well.

Figure 8 Buildings with complex 3D shape

Figure 7 Time-average streamline around building20

coastline. In steep mountain terrain close to the upper right
portion of the coastline, tsunami inundation is obstructed. Somestrong reflected wave and local concentration of high wave
height distribution in narrow steep terrain and inhibitory effectof buildings and structures are confirmed. It also can be seen
how the mountain and building prevent the tsunami, or how thetsunami invasion proceeds along open spaces in detail (Fig.11, 12).
As a result, the 3D tsunami simulation using VOF
method is found as an ef ficient method that can reproduce
the behavior of the tsunami under considering the effects ofsuch building blocks and wide urban area with complex terrain,
nowadays.

Figure 9 An example of a computational mesh

3.3 A large scale 3D tsunami simulation using VOF method [6]Figure 10 Computational domain and tsunami run-up
in a wide urban areamesh

Two-dimensional shallow water model is known as a traditional
tsunami analysis that can be used to predict the inundation
depth and invasion of tsunami in a wide area. However, thismodel could not be used directly in prediction of tsunami force
a c t i n g o n b u i l d i n g o r e v a l u at i o n o f t h e e f f e c t s o f l o c a lsurrounding environment at the level of building in detail.
In the other hand, the gas-liquid t wo -phase f low
simulation using such as VOF method [7 ] can reproduce thewave breaking process well, and has also been applied in somesloshing damper applications of building design. Becauseof the 1000 times densit y changing from the gas to liquid
at the interface, the simulation is known as one of the mostchallenging theme in CFD and it takes a lot of computational
time. Therefore, until now the simulation is impossible to be

Figure 11 Tsunami run-up in complex terrain

used in evaluation of tsunami impacting on a construction ina wide area. Taking into account the performance of computer
which was discussed in the previous chapter, a large-scale
simulation using VOF method is carried including a wide urban
area with complex terrain of 6km*6km. Total number of mesh is
approximately 50million. In order to simulate the 7 minutes timein actual phenomena, it takes 8days in the practical use of using
192CPU cores of the supercomputer TSUBAME 2.0.
As shown in Fig.10, the tsunami is reaching into the

Figure 12 Tsunami run-up in vicinity

21

A Large-scale Simulation on CFD
in Construction Industry

Summary

4

From a practical standpoint, this study has shown that larges c al e simulatio n in massi vel y p ar all el co mp u t atio n is an
important factor to maximize the performance of computeron the CFD applications of the construc tion industr y. The
large-scale simulation is also found as an useful method in the
improvement of calculation accuracy using high resolution
computational grid, in development a more elaborate numericalmodel reproducing a much more flexible shape of an object,and in a wide-area simulation of such 3D tsunami analysis using
VOF method. These new calculation methods are now being
established with the development of computers. Therefore,
developing new evaluation methods due to the new calculation
methods will be an urgent issue, nowadays.
Acknowledgements
The numerical simulations were carried out on the TSUBAME2.0
supercomputer in the Tokyo Institute of Technology supportedby the MEXT Open Advanced Research Facilities Initiative.
References[1]

Top500: http://www.top500.org/[2]

OpenFOAM: http://www.openfoam.com/[3]

P.V. Phuc: Computational Fluid Dynamics using An Open
Source, Journal of Japan Society of Fluid Mechanics, Vol.31,
pp.277-284 (2012)[4]

P.V. Phuc, T. Nozu, H. Kikuchi, K . Hibi: A Large Parallel
Co m p u t at i o n a n d G PU o n CFD, Pro c . o f t h e H i g h
Performance Computing Symposium, pp.75 (2011)[5] P.V. Phuc, T. Nozu, H. Kikuchi, K. Hibi, Y. Tamura: Local Wind
Pressure Distributions on Building with a Setback obtainedby Large Eddy Simulations, Proc. of the 22nd National
Symposium on Wind Engineering, Vol.22, pp.359-364(2012)[6]

P.V. P h u c , M . H a s e b e , I . Ta k a s h i : A St u d y o n T h r e e Dimensional Tsunami Analysis using VOF Method, Annual
journal of coastal engineering, JSCE, Vol. 68, pp.71-75 (2012)[7] C．W．Hir t，et al.: SOLA-VOF: A solution algorithm fortransient fluid flow with multiple free boundaries , Los
Alamos Scientific Laboratory，pp.1-32 (1980)

22

Estimation of strong ground motion anddamage of indoor objects
during giant earthquake
Toshiyuki Masatsuki** Kozo Keikaku Engineering Inc.

Two programs have been por ted to TSUBAME 2 . 0 in order to estimate the damage due to giant ear thquake.
One is the program for seismic wave propagation analysis and the other is for the simulation of the behavior
of indoor objects. The strong ground motions due to the Tokai, Tonankai, and Nankai Ear thquake have been
simulated using a number of GPU. The computation time with GPU was about 10 times faster than not using
GPU. The behavior of indoor objects in a warehouse has also been simulated. The computation time was about
two times as fast as before.

Introduction

1

evaluate the strong ground motions. For now, a simple method
is mainly usedex.[1] due to the limitations of time and cost. Seismic
wave propagation analysis is also one of the evaluation methods.
This method enables to simulate the more realistic seismic
ground motions, especially long period ground motions, than

By the 2011 off the Pacific coast of Tohoku Earthquake, a wide area

simple methods. However, this method requires large amount ofof northeastern Japan was devastated. It is also predicted that acalculation and computation time. Thus, it is not commonly usedwide area of Japan will suffer a great deal of damage in case theexcept important structures. But, the guideline of Architectural

Tokai, Tonankai, and Nankai Earthquake occurs. So, it is important

Institute of Japan[2] recommends using more advanced methods

to estimate the earthquake damages and taking measures against

such as seismic wave propagation analysis, so advanced methods

the earthquake damages has become a critical issue.

may be used more often in the future.

Here we introduce a current situation and problems

Nex t, we show the current situation about indoor

concerning ear thquake damage estimation from the point

damage estimation. There is no detail guideline for overturningof safety structural building design. And then, we show theor dropping objects such as furniture unlike structural building

result of industry trial use project of TSUBAME2.0 in FY2011. In

design. In the past years, especially after the 2011 off the Pacificthe project, we have improved the practical use of the groundcoast of Tohoku Earthquake, we have often noticed the danger

motion estimation and the evaluation of indoor damage due toand the damage of indoor objects. Many investigations havea giant earthquake.

been reported ex.[3] and the various issues about indoor damage

At first, two programs have been ported to TSUBAME2.0

are discussed in some committees. It is possible to estimate

in the trial use project. One is the program for seismic wave

realistic behavior of indoor objects using rigid body simulation

propagation analysis and the other is for the simulation of theconsidering detail information. But, such simulation is rarely usedbehavior of indoor objects. Second, the strong ground motions

because it requires a great deal of time for accurate simulation.

due to the Tokai, Tonankai, and Nankai Earthquake have been

If the computation time is drastically shortened, indoor damage

simulated and the behavior of indoor objects in a warehouse

simulation would be used more often in order to discuss thehas also been simulated.

countermeasures against indoor damage.
In Kozo Keikaku Engineering Inc., we have star ted

Current situation of earthquake
damage estimation2

the project to improve the practical use of the seismic wave
propagation analysis in FY2009. In the first year, we have ported theprogram for seismic wave propagation analysis to Earth Simulator 2
using trial use project by the MEXT. In the second year, FY2010, we
have optimized the program and added new functions. We haveFirst of all, we introduce a current situation and problems about

succeeded to shorten the computation time drastically. In addition,

earthquake damage estimation in terms of safety structural

we have tried to calculate the shorter period ground motions tobuilding design.

do more detail simulation. In FY2011, we have ported the program

Seismic ground motion intensity is apparently critical

to TSUBAME2.0 in order to calculate with GPU, and the simulation

factor for the safety of building. There are various approaches toprogram for the behavior of indoor objects have also been ported.

23

Estimation of strong ground motion anddamage of indoor objects during giant earthquake

Wave propagation analysis

3compiler to use CUDA in Fortran. As a result, the computational
speed with GPU was almost 10 times faster than using only CPU
with -fastsse option (3.1 billion mesh, 100 nodes and 300 GPUs).
Fig. 1. is the weak scaling performance of this program. This
program scales well with the increasing number of GPUs. The

Here we describe the result of the trial use project of TSUBAME2.0

bottle neck of this program is MPI communication time, so it isin FY2011.

required to decrease communication data in order to make thisFirst of all, we describe the overview of the program

program faster. If GPU global memory size becomes larger, it

for seismic wave propagation analysis ported to TSUBAME2.0.

is able to proceed with same scale simulation with less parallel

This program calculates seismic wave propagation by 3D FDM.

number, thus the computation time may become shorter. The

The program had already been able to calculate using a number

largest scale simulation of this project was 12 billion subdomainof nodes with MPI. So, we mainly updated this program tosize with 810GPUs using 270nodes. And its computation time wascalculate with GPU by each process by reference to Nakamura

4 hours for 77600 steps (duration time is about 232 seconds).

et al (2011)[4], who simulate seismic wave propagation using same

Nex t, we describe the seismic wave propagation

method with GPU. The program was compiled by PGI Fortransimulation of the Tokai, Tonankai, and Nankai Earthquake. Fig.

Figure 1

Weak scaling of the program
for wave propagation analysis

Low Pass Filter : Period T >= 2.0sec

Figure 3

Velocity Time History in Nagoya

Figure 2

24

Simulation area and the seismic
source of the Tokai, Tonankai,and Nankai Earthquake2. shows the simulation area and the seismic source of thissurface in Nagoya. The ground shaking started at 50 seconds andearthquake. The seismic source parameters were set by reference

continues to around 300 seconds. Fig. 4. shows the snapshot of

to Central Disaster Prevention Council (2003)[5] and Yamamoto,seismic wave propagation. The contour color shows the velocity on[6]

Yoshimura et al (2012) . The subdomain size is 3.1 billion

the surface ground. The 75 seconds after the fault rapture started,(550x750x60 km), 0.2x0.2x0.2 km mesh size. The delta time is 0.006

the seismic wave arrives to the Osaka Plain and the Nobi Plain, andseconds and the number of the simulation step is 50000 (duration

the ground shaking in the basin area continues for a long time. This

time is 300 seconds). The computable period is more than 2

is the same as observed phenomenon during an earthquake. The

seconds. Fig. 3. shows the velocity time history on the groundvalidation of the accuracy of this simulation is the future work.

Figure 4 Snapshots of the seismic wave
propagation analysis(Contour color: Surface velocity)25

Estimation of strong ground motion anddamage of indoor objects during giant earthquake

Indoor damage estimation

4

Now, we show the overview of the 3D rigid body dynamics
simulation program for indoor damage estimation. This program
can simulate the behavior of falling objects from a rack during
an earthquake such as Fig. 5. Analytical Method [7],[8],[9] , which
treats contacts as constraints, is used in this program to calculatecontact forces between colliding objects.
The brief flow of this calculation is as follows;1.

Detect collisions of rigid bodies

2.

Update the constrains due to the contacts based on the

Figure 5

Simulation of falling objects from racks
during an earthquake

result of the collision detection
3.

Calculate constraint forces including contact forces suchas a normal or friction force

4.

Update velocities and positions of rigid bodies based onequation of motion

First step porting this program to TSUBAME2.0, the program was
modified to handle large scale data. Before the modification,
if the number of objects were over about 1000, the simulation
was stuck. After the modification, it is possible to proceed thenumber of objects is even more than 10000. We used Intel TBB
for the parallelization and enables to calculate with multi-GPU.
The process of collision detection were mainly acceralated
because this process spends a big part of computation time. As
the result, the computational speed of collision detection has
become 3 to 4 times faster than before. In consequence, theFigure 6

Simulation model of indoor objects

whole computation time is about 2 times than before.
Fig. 6. shows the simulation model of indoor objects inwarehouse. There are one meter square pallets showing with redline in Fig. 6., and rigid body assumed 48 jute bags are showed

Conclusion

with white lines on the square pallets. 48 jute bags are one

5

block and they are piled up with various directions. The whole
number of jute bags was about 1400. We have confirmed theWe have ported two programs to TSUBAME2.0. One is the program

model include over 10000 rigid bodies was possible to calculate.

for seismic wave propagation analysis and the other is the program

However it took 2 week s to calculate only several tens offor simulating indoor damage due to giant earthquakes.Then,

seconds, and it is noticed that the calculation becomes unstable

we have shortened the computation time and constructed thewhich has enormous number of objects. Therefore we could notenvironment to do the simulations easier than before. We will

do the simulation which has over 10000 objects.

validate the results of the simulations in the future.

Fig. 7. shows the result of the simulation. The wave
shown in Fig. 3. is used as an input floor response. From thisAcknowledgements

figure, it shows the jute bags fell down and one of the cargos is

26collapsed and flopped to next cargo around 76 seconds passed,

The numerical calculations of this work were carried out on thewhich is the time for strong shaking. The result of this simulation

TSUBAME2.0 supercomputer in the Tokyo Institute of Technology

should be validated in the future.

supported by the MEXT Open Advanced Research Facilities Initiative.

Figure 7 Snapshots of the simulation of indoor
objects in ware house27

Estimation of strong ground motion anddamage of indoor objects during giant earthquake

References[1]

S. Midorikawa and H. Kobayashi: On estimation of strong
earthquake motions with regard to fault rupture, Transactionsof the Architectural Institute of Japan, Vol. 282, pp. 71-81 (1979)[2]

Architec tural Institute of Japan: Generation guide forseismic input motions based on the recent advancement ofground motion studies (2009)[3]

Tokyo Fire Department: The 2011 Tohoku Earthquake Survey(2012)[4]

T. Nakamura, H. Takenaka, T. Okamoto and Y. Kaneda: FDM
simulation of seismic wave propagation for an aftershockof the 2009 Suruga Bay earthquake: Effects of oceanbottom
topography and sewater layer, Bull. Seism. Soc. Am., Vol. 102,
No. 6, pp.2420-2435, (2012)[5]

Central Disaster Prevention Council: “special investigation
council concerning Tounankai and Nankai earthquakes, etc.”(No. 10), (2003).

[6]

Y. Yamamoto and C. Yoshimura: Long-period ground motion
simulation of Tokai-Tonankai-Nankai coupled earthquake
based on large-scale 3D FEM, Journal of structural andconstruction engineering. Transactions of AIJ, Vol. 7, No. 677,
p.1055-1064, (2012)[7] D. Baraﬀ: Analytical Methods for Dynamic Simulation of
Non-Penetrating Rigid Bodies, Computer Graphics Proceedings(SIGGRAPH 89) Vol. 23, pp. 223-232, 1989.
[8]

Kenny Erleben: Multibody Dynamics Animation Lecture 12,
http://www.diku.dk/OLD/undervisning/2005f/101/lecture12.
pdf .

[9]

S. Hasegawa, Y. Tazaki, H. Mitake et al.: “Springhead2”,
http://springhead.info/wiki/index.php .

● TSUBAME e-Science Journal No.8
Published 12/03/2013 by GSIC, Tokyo Institute of Technology ©
ISSN 2185-6028
Design & Layout: Kick and Punch
Editor: TSUBAME e-Science Journal - Editorial room
Takayuki AOKI, Thirapong PIPATPONGSA,
Toshio WATANABE, Atsushi SASAKI, Eri Nakagawa
Address: 2-12-1-E2-1 O-okayama, Meguro-ku, Tokyo 152-8550
Tel:+81-3-5734-2087 Fax: +81-3-5734-3198
E-mail: tsubame_j@sim.gsic.titech.ac.jp
URL:
http://www.gsic.titech.ac.jp/

28

vol.

International Research Collaboration
The high performance of supercomputer TSUBAME has been extended to theinternational arena. We promote international research collaborations using
TSUBAME between researchers of Tokyo Institute of Technology and overseas
research institutions as well as research groups worldwide.

Recent research collaborations using TSUBAME
1. Simulation of Tsunamis Generated by Earthquakes using Parallel
Computing Technique
2. Numerical Simulation of Energy Conversion with MHD Plasma-fluid Flow
3. GPU computing for Computational Fluid Dynamics

Application Guidance
Candidates to initiate research collaborations are expected to conclude
MOU (Memorandum of Understanding) with the partner organizations/departments. Committee reviews the “Agreement for Collaboration” for joint
research to ensure that the proposed research meet academic qualifications
and contributions to international society. Overseas users must observe
rules and regulations on using TSUBAME. User fees are paid by Tokyo Tech’s
researcher as part of research collaboration. The results of joint research are
expected to be released for academic publication.

Inquiry
Please see the following website for more details.
http://www.gsic.titech.ac.jp/en/InternationalCollaboration

8