A Novel Massively Parallel Computing Method for Fluid-Structure
Coupling Analysis
Kazuyasu Sugiyama1,2 , Satoshi Ii3 , Shu Takagi1,2 , and Yoichiro Matsumoto1
1

School of Engineering, The University of Tokyo, JapanRIKEN Computational Science Research Program, Japan3
Department of Mechanical Science and Bioengineering, Osaka University, Japan
2

1

Introduction
Numerical simulations of Fluid-Structure Interaction(FSI) problems have been expected to aﬀord much insight into biological and biomedical systems along withrecent advances in high-performance computing. Unlike
conventionally-used Lagrangian methods with a mesh generation/reconstruction procedure, the authors have developed an Eulerian method using ﬁxed regular grids [1]–[6]. It facilitates simulations of systems involving a large
number of dispersed bodies, and avoids a breakdown in alarge deformation owing to the absence of the mesh distortion problem. In view of computational eﬃciency, themonolithic formulation readily gets great performance outof SIMD processing and makes a workload on each compute core equivalent. It allows us to utilize eﬃcient algorithms cultivated for incompressible ﬂuid ﬂow problems.
Nevertheless, an iterative procedure for solving the pressure Poisson equation makes it diﬃcult to exploit the system performance on a recent scalar supercomputer because
the performance of the solver is limited by the memory
bandwidth and repetitive communications rather than theﬂoating-point operation speed. Following a revived Artiﬁcial Compressibility Method (ACM) [7], the authors havebeen developing a novel algorithm to drastically enhancethe eﬃciency and scalability. We will demonstrate its validity and usefulness for the Eulerian FSI simulation.
2 Numerical methods
2.1 On Eulerian description
As in many analyses for biological systems, the ﬂuid andsolid phases are assumed to be incompressible. The governing equations consist of mass and momentum conservations. The basic equations are discretized on a ﬁxed Cartesian mesh in a ﬁnite-diﬀerence manner. The hyperelastic
dynamics is incorporated into a standard incompressible
ﬂow solver to realize the ﬂuid-structure coupling analysis.
Here, we shall focus on an evaluation of the stress tensor,
which is a key issue in the Eulerian formulation. For moredetailed description, see [2].
In the Lagrangian method, the ﬂuid and solid phases aredistinguished by the boundary ﬁtted mesh, and the level ofsolid deformation is quantiﬁed by the relative change in theadjacent material points. By contrast, in the present Eulerian method [2], to represent the kinematics of structure,
the solid volume fraction φs and the left Cauchy-Green
deformation tensor B are deﬁned on each grid point, andtemporally updated by solving transport equations∂t φs + v · ∇φs = 0,

(1)
L∂t B + v · ∇B = L · B + B · L ,(2)where v denotes the velocity vector, and L(= ∇v T ) thevelocity gradient tensor. The Cauchy stress tensor σ iswritten in a ﬂuid-solid mixture form. For example, a system consisting of Newtonian ﬂuid and neo-Hookean material obeys
σ = − P I + {(1 − φs )μf + φs μs }(L + LT )+ φs G(B − tr(B)I/3),(3)

where P denotes the pressure, μ the viscosity, and G themodulus of transverse elasticity.
2.2

Artificial compressibility method with adaptive parameters
The pressure deviation p from the driving pressure gradient to pump the system is given to satisfy p = 0 (here,... stands for the volume average over the whole computational domain). Following the SMAC algorithm, weintroduce the unprojection and projection steps to update
the variables. The predicted velocity v ∗ is supposed to begiven after the unprojection step. Following [7], we write
the incremental pressure δp(= pN+1 − pN ) as
δp = −(CD∗ + ΓpN )/(Δt),(4)

where D(= ρ∇·v) denotes the divergence of the mass ﬂux.
At the projection step, the velocity is updated asv N+1 = v ∗ − (Δt)ρ−1 ∇(δp).

(5)

Although the solenoidal condition (∇ · v = 0) is not perfectly satisﬁed in the ACM, it is approximated as exactlyas possible. From (4) and (5), we obtain
(DN+1 )2  = (D∗ )2  + 2D∗ ∇2 D∗ C + (∇2 D∗ )2 C 2+ (∇2 pN )2 Γ2 + 2D∗ ∇2 pN Γ + 2(∇2 D∗ )(∇2 pN )CΓ.
(6)
Under the condition of Γ ≥ 0, the parameters C and Γ
are dynamically determined to minimize (DN+1 )2 . Ina MPI parallelization, the advocated algorithm involves
no iterative procedures, and thus considerably reduces theamount of the memory access and adjacent communication. When the second-order ﬁnite diﬀerence is applied,the√ pseudo Mach number was found to have a upper limit
2 3η (here, η is the CFL number based on the maximum
advection speed), guaranteeing the computed velocity ﬁeld
to be nearly incompressible as long as η is suﬃciently
smaller than unity.

Results and discussion
To demonstrate the validity of the method, simulated
results of a channel ﬂow involving 16 discoid biconcave
particles are exempliﬁed. For comparison, two kinds ofsimulations are performed: one is based on the SMAC algorithm, in which the Fast Fourier Transform (FFT) andthe Tri-Diagonal Matrix Algorithm (TDMA) are applied
to exactly solving the pressure Poisson equation, while theother is based on the present ACM. The particle positionand orientation are shown in Fig. 1. The particles deform and translate in the downstream. As the time goeson, they rotate and tend to be more mixed. There are nosigniﬁcant discrepancies in the particle position and shape
between the results based on the SMAC algorithm andpresent ACM. To further check whether an energy conservation is correctly captured, we report a budget of theoverall kinetic-energy transport in Fig. 2. The chained
curve in Fig. 2 corresponds to the summation of all thecontributions, which should be zero. Its absolute value ismuch smaller than the variation of the contributions of theindividual terms. Therefore, the system is well conserved
during the simulation in view of the energy balance, andthe energy exchange between the ﬂuid and solid phases viathe solid deformation is reasonably guaranteed.

Fig. 1.

The snapshots of hyperelastic particles in a threedimensional Poiseuille ﬂow at t = 6.4ms. The computationalextent is Lx × Ly × Lz = 21.12μm × 21.12μm × 21.12μm, andthe number of grid points is Mx × My × Mz = 128 × 128 × 128.
The material properties and the driving pressure gradient are
scaled using the density ρ, and set to μ/ρ = 1 (μm)2 /μs,G/ρ = 2 × 10−2 (μm)2 /(μs)2 and −ΔP/(ρLx ) = 2 × 10−4μm/(μs)2 . The left and right panels show the results based on
the SMAC method using FFT-TDMA and those on the presentACM, respectively. The colors on the walls indicate the magnitude of the shear stress.

Budget [kPa/s]

1.5
1.0
0.5
0.0

I−εs
−εf
−dK/dt

Total

40

30
01010

1210

10

3410

Number of nodes NFig. 3. The ratio of the sustained performance to the peak
performance of the K computer as a function of the number
N of compute nodes. The number of grid points per node is512 × 128 × 128.

4Conclusion
A novel numerical method for ﬂuid-structure interaction systems has been developed. It was applied to thethree-dimensional simulation of the channel ﬂow including hyperelastic particles. Its validity and also high eﬃciency and scalability in massively parallel computing weredemonstrated.
Acknowledgment
This research was supported by Research and Development of the Next-Generation Integrated Simulation ofLiving Matter, a part of the Development and Use of theNext-Generation Supercomputer Project of the Ministryof Education, Culture, Sports, Science and Technology.
Part of the results was obtained by early access to the K
computer at the RIKEN AICS.
References[1] Sugiyama, K., Ii, S., Takeuchi, S., Takagi, S. and Matsumoto, Y., Full Eulerian simulations of biconcave neoHookean particles in a Poiseuille ﬂow, Comput. Mech., 46,147–157, 2010.

[3] Ii, S., Sugiyama, K., Takeuchi, S., Takagi, S. and Matsumoto, Y., An implicit full Eulerian method for the ﬂuidstructure interaction problem, Int. J. Numer. Meth. Fluids, 65, 150–165, 2011.

-1.00.250[2] Sugiyama, K., Ii, S., Takeuchi, S., Takagi, S. and Matsumoto, Y., A full Eulerian ﬁnite diﬀerence approachfor solving ﬂuid-structure coupling problems, J. Comput.
Phys., 230, 596–627, 2011.

-0.5
-1.5
0.0

cellent scalability is found therein. Notably, the performance ratio is 46.6%, 45.6% and 43.7% for N = 1, N = 2
and N = 12, 288, respectively, revealing unprecedentedly
high eﬃciency among numerical simulations of incompressible ﬂuid/structure dynamics performed on state-of-theart scalar-type supercomputers.

flops/Peak (%)

30.4

0.6

0.8

1.0Time [ms]Fig. 2.

The budget of the kinetic-energy transport in thePoiseuille ﬂow containing 16 particles using the present ACM.
Therein, dK/dt, I, εs and εf denote the kinetic-energy rate, theenergy input rate, the strain energy rate, and the kinetic-energy
dissipation rate, respectively.

To demonstrate the applicability of the novel formulation and algorithm to massively parallel computing, theFSI simulations in the channel are also performed for various number N of compute nodes on the K computer. Fixing the problem size per node, we conduct weak scaling
tests, and report the performance in Fig. 3. An ex-[4] Takagi, S., Sugiyama, K., Ii, S. and Matsumoto, Y., A
review of full Eulerian methods for ﬂuid structure interaction problems, J. Appl. Mech., 79, 010911, 2012.
[5] Ii, S., Gong, X., Sugiyama, K., Wu, J., Huang, H. and Takagi, S., A full Eulerian ﬂuid-membrane coupling methodwith a smoothed volume-of-ﬂuid approach, Comm. Comput. Phys., 12, 544–576, 2012.
[6] Ii, S., Sugiyama, K., Takagi, S. and Matsumoto, Y., A
computational blood ﬂow analysis in a capillary vessel including multiple red blood cells and platelets, J. Biomech.
Sci. Eng., 7, 72–83, 2012.
[7] Ohwada, T. and Asinari, P., Artiﬁcial compressibility
method: asymptotic numerical method for incompressible
Navier-Stokes equations, J. Comput. Phys., 229, 1698–1723, 2010.