RIKEN R-CCS Annual Report FY2019

Chapter 13

Complex Phenomena Unified
Simulation Research Team
13.1

Members

Makoto Tsubokura (Team Leader)
Keiji Onishi (Researcher)
Rahul Bale (Researcher)
Koji Nishiguchi (Postdoctoral Researcher)
Kazuto Ando (Technical Staff)
Hsueh-Jui Lu (Student Trainee)
Mehdi Badri (Student Trainee)
Chung-Gang Li (Visiting Researcher)
Ryoichi Kurose (Visiting Researcher)
Tetsuro Tamura (Visiting Researcher)
Shigeru Okazawa (Visiting Researcher)
Huilai Zhang (Visiting Researcher)
Leif Niclas Jansson (Visiting Researcher)

13.2

Overview of Research Activities

The objective of our research team is to propose a unified simulation method of solving multiple partial differential equations by developing common fundamental techniques such as the effective algorithms of multi-scale
phenomena or the simulation modeling for effective utilization of the massively parallel computer architecture.
The target of the unified simulation is supposed to be complex and combined phenomena observed in manufacturing processes in industrial cycles. Our final goal is to contribute to enhance Japanese technologicalcapabilities and industrial process innovation through the high-performance computing simulation.
Most of the complex flow phenomena observed in manufacturing processes are relating to or coupled withother physical or chemical phenomenon such as turbulence diffusion, structure deformation, heat transfer, electromagnetic field or chemical reaction. While computer simulations are rapidly spreading in industry as usefulengineering tools, their limitations to such coupled phenomena have come to realize recently. This is becauseof the fact that each simulation method has been optimized to a specific phenomenon and once two or moresolvers of different phenomena are coupled for such a complicated target, its computational performance is seriously degraded. This is especially true when we utilize a high-performance computer such as Fugaku. In such
125

126

CHAPTER 13. COMPLEX PHENOMENA UNIFIED SIMULATION RESEARCH TEAM

a situation, in addition to the fundamental difficulty of treating different time or spatial scales, interpolationof physical quantities like pressure or velocity at the interface of two different phenomena requires additional
computer costs and communications among processor cores. Different mesh topology and hence data structures
among each simulation and treatment of different time or spatial scales also deteriorate single processor performance. We understand that one of the keys to solve these problems is to adopt unified structured mesh anddata structure among multiple simulations for coupled phenomena. As a candidate of unified data structure forcomplicated and coupled phenomena, we focused on the building-cube method (BCM) proposed by Nakahashi[1].
In summary, the overview of Research Activities in FY2019 was as followings.
In the research and development of unified simulation framework ’CUBE’, an empirical analysis based on
the actual vehicle aerodynamic problem has been conducted. The aerodynamic drag difference was predicted
for the 11 specifications of the small vehicle model, and the results showed a generally good agreement withexperiment. However, there were some cases where the size was reversed. The 16 specifications of SUV typeautomobile models was analyzed, including the rotation conditions of tires with tread patterns. Although it
showed good agreement with the experiment in general, there were cases where the tendency could not be
reproduced. In the future, we plan to continue working to improve accuracy by enhancing practical know-how
together with our partners.
A study of the turbulence model was conducted that is important in the empirical vehicle aerodynamic
analysis, the CSM turbulence model based on the coherent structure of the flow has been implemented. Basedon the results of the turbulence energy decay spectrum of the uniform isotropic turbulence, it was confirmedthat accurate modeling was done. This model has already begun to be applied to practical vehicle aerodynamic
analysis. We are planning to proceed further research to improve its reliability.
As a prediction verification of narrow band noise generated in automobiles, the core technologies of the
CUBE compressible unified solution method has been adjusted to create an improved results. A large-scale
analysis of the real vehicle model has been conducted and successfuly reproduced the acoustic vibration basedon acoustic feedback. In addition, the grid resolution required to capture the correct peak was found fromthe basic verification study, and the required calculation scale has been estimated for the future study. In the
future, we plan to realize analysis with this resolution in “Fugaku” and lead to the creation of results.
The fluid-structure unified analysis method with marker particles has been modified to deal with compressible
solid analysis and verification analyses were conducted. The solid stress distribution in the tensile analysis andthe shear analysis of the circular plate was compared and verified with the reference solution by the finite
element method. The result that the ratio of the cell size to the diameter of the circular hole was 0.78% or lessshowed good agreement with the reference solution. It was also confirmed that the numerical solution by thismethod converged to the reference solution by making the cell size finer.
In the PIC (Particle in Cell) based Euler structural analysis method, the algorithm was applied to the rigidity
analysis of the vehicle white-body by increasing the calculation scale. An analysis of about 210 million cells
was performed on the real vehicle model, and the Mises stress distribution at each time obtained by LS-DYNA
was compared. Qualitatively valid results were obtained, and it was confirmed that the calculation turn-aroundtime could be significantly shortened by using large-scale computational resources. In the future, based on the
results obtained in this research, we will strive to create further results through the activities of the automobile
consortium.
As a deployment of the unified coupled analysis method to practical problems, the dolphin-kick swimming
simulation has been conducted as the field of biomechanics and realized a moving boundary fluid simulation
with smooth motion. By focusing on the width of the legs and the direction of the toes, it was clarified that thedirection of the toes affects the dolphin-kick movement rather than the width of the legs. In the future, we willproceed with analysis of other athletes considering the difference due to body shape, analysis with increased
resolution to improve visualization of vortices. We will continue to cooperate with related research institutes
toward the 2020 Tokyo Olympic Games to create results.
As an application example in the field of biomechanics, an aerodynamic simulation considering the dynamic
attitude change of ski jumping has been conducted, and realized numerical fluid analysis for a series of movementsfrom takeoff movement to initial flight movement. This is probably the world’s first achievement. In the future,
we plan to continue research in order to analyze different jumpers and to solve the problem of clearly explaining
the causal relationship between flight distance, aerodynamic force, and flow structure.
As an application example of the unified compressible solution, a sibilant fricatives pronunciation simulationof human mouse has been conducted. As an oral cavity model, simulations were performed on both a simplified
flow channel and a realistic human model. The results of the sound source spectrum showed good agreement withthe benchmarks in the references. In the future, we plan to clarify the effect of complex shapes on vocalization

13.2. OVERVIEW OF RESEARCH ACTIVITIES

127

and proceed with research to help support denture design.
A performance evaluation and tuning of the incompressible solver has been conducted on K-computer. The
performance improvement of the entire program was obtained as 2.9 times faster at 8,192 nodes than thestandard numerical scheme (Multigrid + CG method) used in the incompressible framework. The speedup issuitable for creating results in “Fugaku”. It has already been confirmed that the unified solution compressible
solver is expected to improve the performance by 12.8 times per node and about 25 times in total in parallel
performance compared to the K-computer by the estimation of the post-K computer performance evaluation
tool. In the future, we plan to evaluate scaling performance further when the number of nodes is increased anddeal with the expected increase in collective communication costs.
As a basic research for constructing a flow solution method using machine learning technology, the reconstruction of a flow mode using a neural network has been tried. When the original flow field was decomposed
into two modes, and reconstructed it again, the flow field could be successfully reproduced with the same accuracy as the reference. We plan to extend this method in the future and evaluate the effectiveness of this method
for 3D problems.
A feasibility study based on some well-known open source code including AMReX with the aim of implementing the adaptive mesh refinement has been conducted. Some basic simulations was conducted in search of
new features suitable for our code, and gained important insights in investigating the low implementation effort
and missing features in the available code. In the future, we plan to continue more detailed investigations andfinally consider implementing it into CUBE.
A multi-objective optimization analysis framework has been developed and a genetic-evolutionary optimization up to the 12th generation with 4-objective functions has been conducted. From the Pareto solution in theobtained aerodynamic drag (Cd) value and the difference in drag (∆Cd) to the wind direction, it was confirmedthat the tendency of the objective function value was similar to the knowledge obtained in the actual design
field. Furthermore, using these results as training-data, a surrogate model construction was tried by machinelearning with 9 design variables as inputs and Cd and ∆Cd as outputs. The Cd value was predicted almost
accurately in the prediction result by machine learning compared with the actual calculation result by CFD,
and it was confirmed that the constructed neural network was appropriate. In the future, we plan to introduce
design variables that reproduce more detailed shapes and improve the degree of freedom, which will lead to thecreation of results in “Fugaku”.
A topology optimization framework for the Lagrangian marker particle method, which is a structural analysis
method for Eulerian analysis, has been developed. The topology optimization results were verified by comparing
them with those of a commercial finite element structural analysis software (LS-DYNA), and the formation of an
optimized shape equivalent to the previous one was confirmed. In addition, topology optimization was performed
on a design domain with 10 million design variables, which covers the white body frame of an automobile, anda frame shape close to the actual shape was obtained. In the future, we will continue to carry out various basic
verification examples and optimization verifications for a large number of load cases that are assumed to be real
problems, as well as to improve the optimization performance through collaboration with industry.
As part of the development of an automobile engine test bed, we performed a fuel-air mixing analysis thatassumed a practical usage scenario. The simulation results of the standard engine shape showed the acceptableconservation of density, pressure and temperature comparing with the theoretical result. The comparison of thedifferent shape of the piston showed the detailed insight of the flow mechanism for the strength of the tumble
and swirl inside the engine. However, the calculation condition is still under the cold flow condition. In next
step, we are planning to model the combustion inside the engine and conduct direct evaluation of combustionenergy efficiency. And also we have started to provide this framework to the industrial user for the evaluationof the usability of the framework. We will continue to improve practicality through the consortium activities.
As a first step of the development of the combustion system test bed, the flow simulation of the jet engine
combustor was performed. It can be inferred that the simulation qualitatively resolves the injector flow, though
there was a qualitative difference at downstream section between the simulation and the PIV data which does
not have depth component. In the next step, we are planning to develop the fuel spraying model and the
combustion model which is specifically designed for the real engine combustors. We will continue to develop
the functionality to reproduce the real state of engine combustion through the consortium activities.
As an issue of the architectural consortium, verification analysis of a test bed based on a unique fluid
structure coupled analysis method was performed. First, the verification analysis for basic tower structure
was performed. The results of wind pressure comparison with the wind tunnel experiments showed a wellagreed profile. It proves the validity of the method that we have developed as a unique fluid-structure strong
coupled unified solution that combines marker particles with the Euler-type structural analysis method. Next,
the numerical stability was examined by performing analysis with various conditions of prismatic structures.

128

CHAPTER 13. COMPLEX PHENOMENA UNIFIED SIMULATION RESEARCH TEAM

The numerical instability occurred near the prisms due to the numerical calculation method of solid strain. It
revealed that it requires introducing a method that can stably calculate the solid strain without depending on
the velocity field. Next, we conducted a simulation of wind resistance analysis of the dome roof of a stadium
structure as an application to practical problems. We could realize the considerable world’s first example ofsolving FSI problems at a high Reynolds number over hundred million by an Eulerian FSI method. The detailed
analysis of the results has not been completed; however we plan to increase the resolution of calculation meshto try to capture smaller flow structures and deformation modes. We will continue to develop the functionality
to produce the practically useful framework for the strong-coupling fluid-structure interaction, and contributeto help users to apply this method into the practical cases through the consortium activities.
By combining aerodynamic analysis with vehicle motion analysis including suspension, steering mechanism,
driver’s operation reaction, etc., an analysis framework that can reproduce the actual driving state has beendeveloped. The maneuver using a hatchback type real vehicle model has been realized and analized. By focusing
on the roll motion and the steering input by the driver, the results consistent with the reports in the humansensory test were obtained. In the future, we plan to carry out higher-resolution analysis and improve versatility
while making improvements toward the creation results of “Fugaku”.
[1] Nakahashi, K., “Building-cube Method for Flow Problems with Broadband Characteristic Length,” In:
Armfield S., Morgan R., and Srinivas K., editors. Computational Fluid Dynamics. (Springer, 2002, 2003),
pp.77-81.

13.3

Research Results and Achievements

13.3.1

Empirical analysis using unified analysis framework CUBE

13.3.1.1

Verification analysis for WLTP in vehicle aerodynamics

The WLTP (Worldwide harmonized Light vehicles Test Procedure) is an international new standard of fuelconsumption test method in the automobile industry, and an evaluation test method for aerodynamic dragusing computational fluid dynamics (CFD) has been established in automobile aerodynamic design. In Europe(and other member countries), a basic agreement on test methods between companies and certification agencies
will be reached in 2019, and full-scale operation is expected to begin in 2020. On the other hand, in Japan,
although the agreement has been obtained between companies, the situation is still in flux as the approval of thecertification agency, that is the Ministry of Land, Infrastructure, Transport and Tourism, has not been obtained.
It is still important for domestic automobile manufacturers to support WLTP, which expands opportunities forusing CFD, they have requested us that we should continue this research activity. However, it may be a phasethat requires more active lobbying to the certification agency. In this phiscal year, with the aim of further
improvment of accuracy, Suzuki Motor Co., Ltd. and Nissan Motor Co., Ltd. provided us with new vehicle
shapes and actual vehicle measurement data in multiple cases, and we conducted research to improve software
practicality. The rate of improvement in accuracy has already slowed down, and more efforts are required forfurther improvement, so continuous efforts are important. In WLTP, it is expected that the aerodynamic dragmeasurement value by CFD is not an absolute value but only a relative value due to the specification difference.
So in this research, the aerodynamic resistance difference ∆CdA[m] due to the specification difference has been
mainly evaluated.
A total of 11 specifications were analyzed for the 2-box type small vehicle model provided by Suzuki Motor
Corporation. Differences in specifications include not only large spoilers, but also consideration of combinationsof optional parts such as mudguards, tow-bars, and door visors. The grid resolution is arranged with a minimumof about 1.5 mm in the necessary places such as near the grill, and adjusting it between about 1.5 to 6 mm
according to the part of the car body. The number of grids is about 100 million, about 15 hours using 657
nodes of K-computer in each case for unsteady calculation of 300,000 time steps (equivalent to about 3 secondsof physical time), and 32 nodes in Oakforest-PACS that took about 43 hours to calculate. In addition, we plan
to reduce the calculation time further in order to improve practicality. Figure 13.1 shows an example of thetime-averaged flow velocity distribution obtained. It can be confirmed that the local change in the flow aroundeach part affects the entire flow field such as the wake of the vehicle body and the flow in the engine room, andaffects the aerodynamic performance.
Figure 13.2 shows the drag difference (∆CdA) due to the specifications of the obtained aerodynamics dragcoefficient. Although the tendency is generally grasped, there are some cases where the tendency is reversed.
Consideration of grid resolution dependency, aerodynamic drag accumulation graph, reproducibility of windtunnel test, knowledge of analysis results by other commercially available software, consideration of deformation

13.3. RESEARCH RESULTS AND ACHIEVEMENTS

129

for soft materials, CAD shape and actual vehicle shape comparison etc., the boundary conditions have been
readjusted by Immersed Boundary and the turbulence model has been improved. Thus, some improvements
have been obtained. However, it cannot be said that it has sufficient accuracy yet, so we plan to continue
working to improve the accuracy in next year.

Figure 13.1: Profile of time-averaged dimensionless velocity magnitude in central section for several specificationsof small vehicle simulation.

Figure 13.2: Predicted drag delta for all specifications of small vehicle simulation.
A total of 16 specifications were analyzed for the SUV type automobile model provided by Nissan Motor Co.,
Ltd. The actual measurement is conducted using a wind tunnel with a moving belt, and the analysis adopted
that the condition involves tire/wheel rotation. Differences in specifications include a variety of aerodynamic
add-on parts such as spoilers, deflectors, baffles, and grill covers. The grid resolution is adjusted to a minimumof about 1.5 mm near the grill and about 1.5 to 6 mm depending on the part of the car body. The numberof analysis grids is about 120 million. Unsteady calculation of 150,000 time steps (equivalent to about 1.5
seconds of physical time) in each case required about 48 hours using Oakforest-PACS 32 nodes. Tire rotation
provides a velocity boundary on the tire/wheel surface. In addition, the tire model is produced by a 3D scan
technology which has the real shape, and the detailed shape including the tread pattern is reproduced. We
succeeded in reproducing the rotation of the tire shape including the tread, but under this condition there
are still problems that some calculations become unstable (the moment when the calculation point switches
inside and outside the tire, etc.), so we will continue to improve. When this is completed, it will be the world’s
first case. Figure 13.3 shows an example of the time-averaged flow velocity distribution obtained. Figure 13.4
shows the difference (∆CdA) due to the specifications of the obtained aerodynamic drag coefficient. Although

130

CHAPTER 13. COMPLEX PHENOMENA UNIFIED SIMULATION RESEARCH TEAM

it shows good agreement with the experiment, there are cases where the tendency cannot be reproduced. In
particular, since it was difficult to reproduce aerodynamic parts that affect the underfloor flow, we are reviewing
the calculation conditions such as the modeling method of the wind tunnel floor surface and consideration ofshape deformation. In this process, we are also preparing a group of post-treatment tools that contribute toflow evaluation, such as measuring surface pressure and creating an accumulation graph of aerodynamic drag.
We will continue to work on improving accuracy in next year.

Figure 13.3: Profile of time-averaged dimensionless velocity magnitude in central section for several specificationsof SUV vehicle simulation.

Figure 13.4: Predicted drag delta for all specifications of SUV vehicle simulation.

13.3.1.2

Narrow-band noise analysis around the engine hood of the vehicle

The narrow-band noise in automobile industries is caused by the quality problems and production issues andgenerally, it takes tremendous effort to solve this kind of problem in the design site. To improve the prediction
accuracy, that the commercial software couldn’t achieve by using the conventional methods, shorten the development period, and reduce the cost, the unified compressible solver installing in CUBE is under development.
In this fiscal year, the knowhow and core technologies of the compressible solver have been summarized and the
production of the results is under preparing.
For the resolution of the mesh, the finer resolution with a minimum width of 0.2 mm is located near the gapof the hood so that the acoustic feedback generated from the engine hood can be accurately captured. The total
number of the grid is about 2 billion and 13,086 nodes in K-computer for about 50 hours are consumed. Figure13.3. RESEARCH RESULTS AND ACHIEVEMENTS

131

13.5 shows the contour of the propagation of the pressure fluctuation and the propagation with a circle shapeemitting from the hood can be observed. This phenomenon is completely different from the usual acoustic
noise caused by the vortex around the vehicle body. This kind of acoustic feedback can only be generated
when the resonance appears. To capture this phenomenon, besides the much finer resolution compared withusual fluid simulation, a numerical scheme with extremely small dispersion and dissipation are also required.
The resolution is determined by the flow velocity characteristics (the flow velocity, the gap of the hood, andthe cavity length in the engine room), hence, it is necessary to perform the grid test to validate the current
framework.

Figure 13.5: The propagation of the Pressure fluctuation.
In the second half of this year, because of the lack of computational resources, the large scale simulation
couldn’t be smoothly performed. Therefore, the basic validation of the cavity aeroacoustics problem has been
conducted first. Figure 13.6 shows an example of the grid-dependent analysis of the frequency response for thesimple cavity problem. From this result, the resolution required to capture the correct peak can be known,
and the required computational resource can be also estimated. It can be seen that this simulation requires a
resolution of around 0.5 mm (for a cavity width of 20 mm). Therefore, in the analysis equivalent to the realvehicle, it can be known that a resolution of 0.05 mm is required for the hood gap of 5 to 6 mm. In summary,
10 times finer than the current resolution is necessary. In the future, we would like to perform such large scale
simulation in Fugaku and it is expected that it will be the largest scale simulation all over the world by thetime to help us to understand the mechanism of the narrow-band noise.

Figure 13.6: The comparison of the Frequency.

13.3.1.3

Research on turbulence models that contribute to large-scale aerodynamic analysis

In high Reynolds number turbulence analysis represented by automobile aerodynamics, it is important to selecta suitable turbulence model for large-scale analysis. The Sub-grid-scale (SGS) model based on Smagorinsky(Smagorinsky, 1963) is widely used in the transient Large Eddy Simulation (LES), but it may be necessary
to calculate the distance from the wall surface to the grid point (Van driest damping term of the standard
Smagorinsky turbulence model). Or it requres operations that filter the model coefficients in a uniform direction
to stabilize the solution (least squares operation of the dynamic Smagorinsky model [1]), and these methods
has a high incidental calculation costs. In order to avoid this problem, in this study, we have mainly used
the so-called “Implicit LES”, which compares the artificial viscosity introduced by the stabilization method of

132

CHAPTER 13. COMPLEX PHENOMENA UNIFIED SIMULATION RESEARCH TEAM

the numerical scheme with the physical viscosity and we have utilized it into practical simulation cases. In
this year, we implemented a coherent structure turbulence model (CSM) based on the dominant vortical flow
structures, which was derived by Kobayashi et al.[2]. In this model, in modeling vortex viscosity, the local
constants of the Smagorinsky SGS model are dynamically determined based on the relationship between thesecond invariant of the velocity gradient tensor and the energy dissipation rate of the vortex. Unlike the dynamic
Smagorinsky model, it does not require filtering and does not require distance calculation, which is advantageous
for large-scale calculation and application to complicated shapes. The SGS model by Smagorinsky is expressed
as follows:¯ 2 |S̄|,νt = C ∆|S̄| =

p
Sij Sij ,¯ = ∆¯1 ∆¯2 ∆¯3 1/3 is the filter width, given by the grid width ∆¯i in the i-th direction. In the CSM,
where ∆
the model coefficient C is determined as follows. The second invariant of the velocity gradient tensor Q isrepresented as:
1(Wij Wij − Sij Sij ) ,
2

p
1
T
W̄ =
∇ū − (∇ū) , |W̄ | = Wij Wij ,2
Q=

where W̄ is the vorticity tensor. This definition is termed the Q-criterion, and it is often used to visualize
vortex structures [3]. Kobayashi determined the model coefficients (in non-rotating flow) independently of theflow properties and Reynolds number, as follows:
FCS = Q/E,
1
E = (Wij Wij + Sij Sij ) ,2
C = C1 |FCS |3/2 , C1 = 1/20.
As a result, an eddy viscosity model is expressed, in which the coefficient is automatically zero in laminar
flow, and is naturally damped in the wall direction. For details, refer to [2].
Figure 13.7 shows the basic verification results for this turbulence model on the isotropic turbulence decayenergy spectrum. The homogeneous incompressible isotropic turbulence was generated in a three-dimensional
cubic box using the random Fourier mode. The side length of the box was L = 0.09 × 2π m, with a numberof grid (cell) of 128 × 128 × 128 in three dimensions. The number of Fourier modes was 5000 and a genericdynamic viscosity ν = 1.0 × 10−5 m2 /s was chosen. The initial turbulence velocity field was generated fromthe energy spectrum table of Comte–Bellot et al. (CBC) [4], based on the method of Saad et al. [5]. A Fourier
series of an arbitrary spatial velocity field u at point x is given by:u(x) = 2MX
m=1

qm cos(κm k̂m · x + ψm )σ̂ m ,where M is the number of modes, qm is the amplitude, κm is the m-th wavenumber, and k̂m ≡ (kx,m , ky,m , kz,m )is a unit direction vector. The procedure
consisted of choosing random k̂m and ψm with qm obtained from thepenergy spectrum such that qm = E(κm )∆κ. The modal direction vector σ̂ m was evaluated enforcing thedivergence-free constraint; it was selected to be orthogonal to k̂m ,k̂m · σ̂ m = 0,

∀ m ∈ {0, 1, ..., M }.

In Fig. 13.7, the input spectrum corresponds to the CBC data (solid lines). The simulated data are shownin filled circles and correspond to non-dimensional times of t∗ = 42, 98, 171, respectively. The Nyquist limit isindicated via the vertical dashed line. The results of the initial spectrum and the decaying spectrum both well
reproduced the results of the direct numerical simulation, confirming that accurate modeling was performed.

133

13.3. RESEARCH RESULTS AND ACHIEVEMENTS

This model has already begun to be applied to practical vehicle aerodynamic analysis. From the results ofinitial studies, it has been confirmed that reasonable results with more suppressed numerical viscosity can beshown in areas that were previously sensitive to flow separation. In the next step, we plan to continue research
to improve reliability by implementing a wide range of verification calculations and accumulating empiricalcalculation experiences.

1.E-03E(k) [m3/s2]

1.E-04

1.E-05t*=42t*=98t*=171
CBC t*=42
CBC t*=98
CBC t*=171

1.E-06

1.E-0710100
k [1/m]

1000

(b)

(a)

Figure 13.7: Isotropic turbulence decay energy spectrum for the CSM turbulence model. (a) The turbulence
structures visualized by the Q-criterion. (b) The black circle represents the initial reference energy E(k) for
each wavenumber k; the red or blue colored circle shows simulation result. The solid line shows CBC data. The
Nyquist limit is indicated via the vertical dashed line.

13.3.1.4

References[1] M. Germano, U. Piomelli, P. Moin, W. H. Cabot, A dynamic subgrid-scale eddy viscosity model, Physics of
Fluids A: Fluid Dynamics 3 (7), (1991), 1760-1765.
[2] H. Kobayashi, The subgrid-scale models based on coherent structures for rotating homogeneous turbulence
and turbulent channel flow, Physics of Fluids 17 (4), (2005), 045104.
[3] M. Tanaka, S. Kida, Characterization of vortex tubes and sheets, Physics of Fluids A: Fluid Dynamics 5 (9),(1993), 2079-2082.
[4] G. Comte-Bellot, S. Corrsin, Simple Eulerian time correlation of full-and narrow-band velocity signals in
grid-generated “isotropic” turbulence, Journal of Fluid Mechanics 48 (2), (1971), 273-337.
[5] T. Saad, D. Cline, R. Stoll, J. C. Sutherland, Scalable Tools for Generating Synthetic Isotropic Turbulencewith Arbitrary Spectra, AIAA Journal 55 (1), (2017), 327-331.

13.3.2

Verification analysis of Eulerian finite volume formulation using Lagrangianmarker particles for compressible solid

13.3.2.1

Development of Eulerian finite volume formulation using Lagrangian marker particles
for compressible solid analysis

In this study, in order to deal with large-scale structural analysis problems that reproduce complex shapesthat are difficult to calculate with conventional methods, the Eulerian fluid-structure unified analysis method
using Lagrangian marker particles has been developed. In this fiscal year, to deal with structural analysis
problems with compressible materials, improvement of the fluid-structure unified analysis method with marker
particles and implementation of an Eulerian finite volume formulation using Lagrangian marker particles forcompressible solid analysis has been conducted. At this time, the conservation of mass and the equation of
motion of compressible materials are expressed by the following equations [1].

ρρJ = ρ0
∂v+ (v · ∇) v = ∇ · σ + ρb
∂t

134

CHAPTER 13. COMPLEX PHENOMENA UNIFIED SIMULATION RESEARCH TEAM

Here, ρ0 is the mass density at initial condition, ρ is the mass density, J is the volume ratio, v is the velocity, σis the Cauchy stress, and b is the body force. In this study, the solid is modeled as a compressible neo-Hookean
solid and the constitutive equation is given by [2]1
σ = J −5/3 G B − (trB) I + K (J − 1) I,
3
where, G is the shear modulus, B is the left Cauchy-Green deformation tensor, I is a second-order unit tensor,
and K is the bulk modulus.
In order to confirm the validity of the proposed method, the tensile analysis and the shear analysis of aflat plate with a circular hole were conducted and the solid stress distribution is compared with the reference
solution obtained by the finite element method. Furthermore, by verifying the stress distribution when themesh resolution is changed, the mesh resolution required to obtain the desired accuracy is clarified. In thefollowing examples, the structure has Young’s modulus and mass density of the steel, and the void region is3
given ρ0 = 1.0kg/m to avoid numerical instability. The height of the flat plate(structure) is 1.0 m and thewidth is 0.5 m, and a circular hole with a diameter of 0.25 m is arranged in the center of the flat plate. The
length of the plate in the direction perpendicular to the paper surface (y-direction) is assumed to be a unit
length, and the plane strain state is assumed by setting the velocity in the y-direction to zero.
13.3.2.2

Verification analysis by the tensile analysis and the shear analysis of a flat plate with acircular hole

The lower end of the flat plate fixes the velocity component in all directions to zero, and the upper end of theflat plate gives a tensile velocity that increases linearly in the time direction in the positive z-axis direction.
No restraint conditions are given to the other end faces. In the proposed method, the above velocity boundary
conditions are given by the marker particles arranged in the region. In this example, the analysis time is up to
1 ms. Therefore, the displacement of the upper end of the flat plate is 1 µm, and this example can be regarded
as a minute deformation problem. In the calculation by the proposed method, the time increment size is set sothat the CFL number for the stress wave velocity is 0.1874 at each spatial resolution. In order to obtain thereference solution, a sufficiently fine finite element mesh is used, and the finite element size of the orthogonallydivided region excluding the circumference of the circular hole is 1/640 m.
Figure 13.8 shows the reference solution at t = 1 ms, the marker particles at each spatial resolution, andtheir Mises stress distribution. In this calculation, the number of marker particles in one cell is fixed at 2 × 2.
The result that the ratio of the cell size to the diameter of the hole is 0.78% or less shows good agreement withthe reference solution, and it can be seen that the stress concentration at the upper part and the lower part of
the hole is particularly reproduced. In addition, in order to quantitatively verify the validity of the numericalsolution, the time history of strain energy is compared with the reference solution. As shown in Figure 13.9, it
can be confirmed that the numerical solution by this method converges to the reference solution by making thecell size h finer. In addition, Figure 13.10 shows the error  of the strain energy with respect to the reference
solution at t = 1 ms and its spatial convergence. It is confirmed that the error converges at the order of 0.025
and that the error is 0.27% at the cell size h = 1/1024 m.
13.3.2.3

Verification analysis by shear analysis of a flat plate with a circular hole

The analysis conditions are the same as the example in the previous section, except that the upper end of theflat plate is given a shear rate that increases linearly in the time direction in the positive x-axis direction. Figure13.11 shows the reference solution at t = 1 ms, the marker particles at each spatial resolution, and their Mises
stress distribution. The result that the ratio of the cell size to the diameter of the hole is 0.78% or less shows agood agreement with the reference solution, and it can be seen that the stress concentration around the hole isparticularly reproduced. From Figure 13.12, it can be confirmed that the numerical solution by the proposed
method converges to the reference solution by making the cell size h finer. Furthermore, Figure 13.10 also showsthe error of strain energy with respect to the reference solution at t = 1 ms and its spatial convergence. It is
confirmed that the error converges at the order of 0.077 and that the error is 0.17% at the cell size h = 1/1024m.
When the ratio of the cell size to the diameter of the circular hole is 1.56% or more, it can be seen thatthe oscillation of stress occurs at the top corners of the flat plate in Figure 13.8 and Figure 13.11. This isbecause the velocity field is discontinuous between the cell to which the velocity boundary condition is given
and the cells around it. Thus, the velocity gradient calculated from the velocity of these cells is overestimated.

135

13.3. RESEARCH RESULTS AND ACHIEVEMENTSa) Reference(FEM)b) dx=1/64 m(dx/D ≃ 6.25%)c) dx=1/128 m(dx/D ≃ 3.13%)d) dx=1/256 m(dx/D ≃ 1.56%)e) dx=1/512 m(dx/D ≃ 0.78%)f) dx=1/1024 m(dx/D ≃ 0.39%)

Figure 13.8: Tensile deformation problem of a flat plate with a circular hole: Distribution of Marker particles
and their Mises stress distribution at t = 1ms.

Figure 13.9: Tensile deformation problem of a flat plate with a circular hole: Time history of strain energy.
In order to avoid this numerical oscillation even when the spatial resolution is low, it is necessary to introduce
a calculation method for solid deformation tensors that does not require velocity gradient calculation and wecontinue to work on development to improve analysis accuracy.
13.3.2.4

Large-scale real vehicle model structure analysis using Euler structure analysis framework

In this research, in order to conduct a large-scale structural analysis that reproduces detailed shapes of realvehicle model that are difficult to calculate with conventional methods, we are developing a unique fluid-structure
unified solution method that combines an Euler-type structural analysis method with marker particles. We
applied the algorithm to the stiffness analysis of the vehicke white-body by increasing the scale of calculation.
3
In this report, steel is assumed as the solid material in all the numerical analysis examples, and ρ = 7850kg/m ,E = 200GPa, µ = 0.3 are given as the material parameters of the solid. Here, E is Young’s modulus and µ3
is Poisson’s ratio. In the region where no solid exists, ρ = 1.0kg/m is set to avoid zero division in numericalcalculation. The shape of the automobile body and the hierarchical Cartesian grid based on the building cube
method are shown in Fig. 13.13. The minimum cell size was 1.86 mm, and the total number of cells was about
210 million. 32 × 32 × 32 cells are evenly spaced in a cube. As shown in Fig. 13.14, the boundary conditions in
the torsional stiffness test are reproduced by setting the velocity for the marker particles shown in red. The twomarker particle groups at the rear of the vehicle body are completely fixed, and the two marker particle groups

136

CHAPTER 13. COMPLEX PHENOMENA UNIFIED SIMULATION RESEARCH TEAM

Figure 13.10: Error of strain energy at t = 1ms with respect to reference solution (finite element method) andits spatial convergence.

a) Reference(FEM)b) dx=1/64 m(dx/D ≃ 6.25%)c) dx=1/128 m(dx/D ≃ 3.13%)d) dx=1/256 m(dx/D ≃ 1.56%)e) dx=1/512 m(dx/D ≃ 0.78%)f) dx=1/1024 m(dx/D ≃ 0.39%)

Figure 13.11: Shear deformation problem of a flat plate with a circular hole: Distribution of Marker Particle
and its Mises Stress Distribution at t = 1ms.
at the front of the vehicle body give a velocity that increases linearly in the positive and negative directions of
the z-axis. In this example, the time increment was set to 6.25 × 10−5 ms, and calculations were performed up
to 4800 steps (t = 3.0 ms). 13,065 nodes of K-computer were used for this calculation. In order to compareand verify the numerical solutions by the proposed method, the reference solution was calculated using thecommercial solid-structure analysis code LS-DYNA. The total number of elements of the finite element mesh in
LS-DYNA was about 660,000, and it is consisted with about 580,000 shell elements, about 10,000 solid elements,
14 beam elements, and about 70,000 rigid body elements.
Figure 13.15 shows the Mises stress distribution at each time obtained by the proposed method and LSDYNA. However, the color bar of Mises stress is normalized. From these figures, it could be confirmed that theMises stress distribution obtained by the proposed method is qualitatively generally valid. Table.13.1 shows theturn-around time of the proposed method and LS-DYNA. The turn-around time could be significantly shortenedin the calculation by the proposed method by using a large-scale computational resource (K-computer 13065
node × 2.4 hours). In the future, based on the results obtained in this research, we will continue to create
further results through the activities of the automobile consortium orginized by RIKEN.
13.3.2.5

References[1] G. A. Holzapfel, “Nonlinear Solid Mechanics: A Continuum Approach for Engineering,” Wiley, 2000.
[2] R. W. Ogden, “Non-linear elastic deformations,” Courier Corporation, 1997.

137

13.3. RESEARCH RESULTS AND ACHIEVEMENTS

Figure 13.12: Shear deformation problem of a flat plate with a circular hole: Time history of strain energy.

Figure 13.14: Geometric boundary condiFigure 13.13: Overview of Euler grid and vetions: (red) velocity boundary on a jig given
hiclewhite body shape (finest size: 1.86mm,
by marker particles.
number of grids 210 million).

Table 13.1: Turn-around time comparison
Grid generation
Calculation
Visualization
Total (Turn-around time)

LS-DYNA
several weeks
0.75 hours
1 hourseveral weeks

Proposed method
10 minutes
2.4 hours
1 hour
3.6 hour138

CHAPTER 13. COMPLEX PHENOMENA UNIFIED SIMULATION RESEARCH TEAM

Figure 13.15: Analysis results by CUBE (left) and LS-DYNA (right) (Mises stress distribution).

13.3. RESEARCH RESULTS AND ACHIEVEMENTS

139

13.3.3

Biomechanics CFD–Development of unified coupled analysis methods into
practical problems

13.3.3.1

Role of feet orientation on thrust generation in dolphin-kick swimming

In swimming competitions using the 50m pool, under water swimming is permitted up to 15m from the starting
point and after turning back all swimming styles other than breaststroke. Generally, the swimming style
employed for underwater swimming is known as dolphin kick. As a result, there is much interest in optimizing
the swimming performance of the dolphin kick swimming.

Figure 13.16: The four feet orientations considered in this study.
Since it is difficult to measure the fluid forces generated unsteady swimming such as dolphin kick swimming
experimentally, various studies using computational fluid dynamics have been conducted. The results of suchstudies suggest that lower limbs and feet account for major contribution to the propulsive forces, and that thejoint angle of the lower body is important. However, most of the prior studies focus on a standard orientationof the legs and feet. Therefore, in this study, we investigate the role of feet orientation on the propulsive forcesin dolphin kick swimming using numerical simulations. For that purpose, the position of the legs such as thewidth of the legs and the direction of the ankles was changed variously using the body shape data of the athlete
measured using the 3D scanner. Experimentally obtained dolphin-kick motion was employed of the simulations.
The movement of the lower body was converted into coordinates and given as motion data. The CAD model ofthe swimmer’s body was coupled with the motion capture data to generate CAD model based motion model.

Figure 13.17: Instantaneous thrust generated by the four feet poses.
Numerical simulation was carried out using this setup for four different feet orientations, the results ofwhich are shown in Fig.13.16. From the results of these simulations we find that there are differences in theflow velocity distribution, the vortex structure generated on the sole of the foot when the leg is kicked down,
and the flow velocity vector. By comparing these results focusing on the width of the legs and the directionof the toes, it was confirmed that the direction of the toes rather than the width of the legs affects the diving

140

CHAPTER 13. COMPLEX PHENOMENA UNIFIED SIMULATION RESEARCH TEAM
Table 13.2: Comparison of averaged thrust produced by the four feet orientations.
Pose
Thrust (N)

A
213

B
220

C
181

D
134

method movement. It was confirmed that collecting the flow velocity vectors generated on the sole of the footmore centrally and posteriorly with the toes facing the center has a large effect on propulsive forces.

Figure 13.18: Visualization of the flow around the feet at the end of the power stroke.
Next, the thrust generated by the swimmer is plotted compared in Fig.13.17. The vertical axis is the forceof the swimmer, and the horizontal axis is the time. For each leg orientation, the red circled area where theswimmer has the most propulsion. From the figure, it can be seen that Pose A and Pose B produce more thruston average compared to the other two feet orientations. The time averaged thrust generated by each of the four
feet orientations are shown in Table. 13.2.
PoseB has the most propulsion power, and PoseA also has a value close to PoseB, but PoseC and PoseD
have relatively low values, and PoseD has the lowest propulsion power. From Fig.13.18, it can be seen that theflow velocity vector is generated in a circle on the sole of the foot.
In summary, we simulated the diving movement by intentionally changing the leg position for one athlete,and compared the effect of the leg position on the propulsive force. The following findings were obtained fromthe simulation results.
1. There was a difference in the flow velocity distribution and the force generated by the swimmer dependingon the position of the legs.
2. The magnitude of the vorticity generated on the sole of the foot when the leg was kicked down greatly
differed depending on the width of the leg opening and the direction of the toe.
3. There was a correlation between the vortex generated on the sole of the foot and the magnitude of theflow velocity vector.
4. The direction of the toes had the greatest effect on propulsion, followed by the width of the legs that had
the greatest effect on propulsion.

13.3. RESEARCH RESULTS AND ACHIEVEMENTS

141

From the above, it is considered that the position of the legs in the diving swimming method can be more
propulsive by the posture in which the legs are closed and the toes are aligned with each other as in Pose B.
In optimizing the diving method, as a future research subject, it is necessary to consider the difference due tobody shape, perform the same analysis for other athletes, and perform more precise analysis with increased
resolution to improve the visualization of vortices.
13.3.3.2

Aerodynamic simulation of the dynamic flight posture in ski jumping

In ski jumping competitions, take-off and flight posture are important movements that are greatly influence theflight distance. These two postures are closely related, and errors in the take off posture will adversely affectthe flight posture and it is usually quite difficult to make adjustments during the flight. On the other hand,a jumper must dynamically adjust the flight posture taking into account the wind condition to prevent theincrease in drag and the decrease in lift as much as possible. Failure to do so can result in the jumper stalling
mid-flight, greatly reducing the jumping distance. Although these two movements are very important, there
has been no study that aerodynamically investigated both movements together. This is because it is difficult tocontinuously reproduce both movements in a wind tunnel experiment, and it is also difficult to reproduce a large
flight posture change in simulations that employ unstructured grid that that is typically used for such analysis.
However, with the advent of local mesh refinement techniques in cartesian meshes and immersed boundary
methods that can faithfully reproduce the dynamic motion of a ski-jump, for example, it is now possible tocarryout such simulations. Therefore, the purpose of this study is to carry out aerodynamic simulation using
the above framework for take-off motion to initial flight motion, and to investigate the characteristics of the
aerodynamic force on the jumper. In this study, in order to reproduce more detailed movements before analysis,
we asked for the cooperation of an active ski jumper and attached a motion sensor to the jumper to measure
posture change data. In addition, a GPS device was attached to the jumper to measure the flight trajectoryand instantaneous speed. Using these data, aerodynamic simulations were performed for movements from 0.3
seconds before takeoff to 0.7 seconds after takeoff.

Figure 13.19: (Left) Schematic of ski jump. (Right) Computational domain setup.
The ski jumping competition is one of the Nordic skiing competitions held in the stadium, and is competedin two elements: distance and beauty of flying type. In ski jumping competitions, a series of movements can
be classified into six categories: run-up, take-off, initial flight, stable flight, landing preparation, and landing(Fig.13.19). Of these movements, two important elements for extending the flight distance are take-off andflight posture (initial flight / stable flight). In the take-off motion, while maintaining the speed close to 90 m /s obtained during the approach run, the posture is changed significantly from the crouched position to the flight
posture within a very short span of time, from 0.2 to 0.3 s. This movement, which involves a large change inposture in a short period of time, is considered to be a very difficult and important movement. In flight, jumper
takes a posture to obtain a large lift and low drag. Although there is no significant change in posture during
flight, it is necessary to fine-tune the posture taking into consideration of the wind conditions. In ski jumping
competitions, once a stall occurs during flight, it is difficult to extend the flight distance, so stall prevention is
also very important.
Table 13.3 shows the start gate, flight distance, Run time, and flight time of each trial measured in this142

CHAPTER 13. COMPLEX PHENOMENA UNIFIED SIMULATION RESEARCH TEAM
Table 13.3: Information obtained in each jump.
First jump
Second Jump
Third Jump

Starting Gate
26
17
22

Flying Distance (m)
112
90
113

Run time (s)
6.1
5.7
5.8

Flight time (s)
4.2
3.2
4.1

study. In the measurement in this study, the largest gate is No. 26 and the smallest gate is No. 17, so there
is a difference of 4.5 m in the maximum approach distance. In this study, the boundary conditions were set asshown in Fig.13.19.
From Fig. 13.20, in the first attempt, there was no significant increase in drag immediately after takeoff(around 0 s at time), which was seen in other attempts, and the drag decreased significantly between 0.3 s and
0.4 s after that. On the other hand, in the 2nd and 3rd trials, the drag increases immediately after the takeoff,then decreases and increases once, and then decreases again. At around 0.7 s time, there was no big difference
between each attempt.

Figure 13.20: Evolution of drag, lift, lateral force and lift-drag ratio.
First, as with drag, lift has a low value before takeoff (before time 0 s) and tends to increase significantly
by taking a flight attitude after takeoff (after time 0 s). However, in the first attempt, the lift was very high attime -0.2 s, and a fall similar to drag was seen at time 0.3 s to 0.4 s. On the other hand, it can be seen that in
the second attempt, unlike the other attempts, the high value was maintained between the times of 0.3 s and
0.5 s. In the third attempt, there was no significant decrease after takeoff, but as with drag, it continued todecrease gradually At 0.7 s, as with drag, there was no significant difference in lift between the three attempts.
The lift to drag ratio is defined by the following equation.
L/D =

Lift
Drag

(13.1)

Regarding the time transition of the lift-drag ratio, for each attempt, the value was close to 0 before thetakeoff (before time 0 s), and after take-off L/D rapidly increases from 0 to 1 within a span of 0.1 s. As for thedifference between the trials, it can be seen that the lift-drag ratio increases significantly only in the first run
at around 0.4 s. A snapshot of the streamlines of the flow at 0.6s after take-off is presented in Fig.13.21.

13.3. RESEARCH RESULTS AND ACHIEVEMENTS

143

Figure 13.21: Streamlines at t = 0.6s for first, second and third jump, respectively.
13.3.3.3

Sibilant fricatives pronunciation simulation using a unified compressible solution

The simulation of human speech organ contributes to the dental application. For example, a well-designed
denture enables patients speak better, especially improvement in pronouncing words containing sibilants orfricatives. However, human phonation process is a tough objective for the simulation because the fluctuationpressure which generates the sound is much smaller compare with the atmospheric pressure. Specifically, theamplitude of sound wave is only 0.01% of standard atmosphere. That means the simulation system is sensitive tothe accuracy of the numerical scheme. Therefore, a numerical framework which is suitable for high-performancecomputing is built for modelling the sibilant fricative production in this study.
The implicit time scheme with immersed boundary method based on a hierarchical structure grid as well as
the modified solution-limited time stepping method for the compressible flow are adopted. Firstly, the acoustic
resonance generated from the flow around the single plate is simulated to validate the numerical scheme andthe result shows that this framework is highly efficient and suitable for the massive parallelization system totremendously save the calculation time. Then, the simulation for a simplified model of the sibilant /s/ isconducted and SPL profiles are in good agreement with the experimental results. Finally, the simulation fora realistic geometry of sibilant /s/ scanned from the human vocal tract is performed to demonstrate that thisframework is capable of making a contribution to the dental application, such as helping designing the dental
prosthesis.
The LUSGS implicit algorithm and the Unified Time Stepping scheme are applied for the time advancement.
The Roe scheme with 5th order MUSCL is used to approximates the convective terms. And the second order
central difference method is adopted to calculate the magnitudes of the viscous terms. The absorbing boundary
condition is used at the outlet for avoiding the computational domain being polluted by the reflecting pressure
waves. The immersed boundary method with the hierarchical structured grid is used to simulate the complex
geometry surface. And the OpenMP and MPI are applied for parallel computation.
The resonance plays an important role in acoustics simulation, since it enhances the amplitude of soundwave and makes characteristics of sound clearly. Hence, the model of single plate is simulated to check thephenomenon of resonance with the benchmark [1] by present framework. In Fig. 13.22, the fluctuation pressure
is indicated. The resonance occurs in a half-wavelength mode along the plate length at U = 44m/s, both of the
wave length and intensity of the sound wave can be observed.
The sibilant fricative is the common consonant in each language produced by forcing air through a narrow
channel inside the oral cavity. Nozaki et al. [2] replicated a realistic model of a subject pronouncing /s/ andused incompressible solver to simulate the flow field. Yoshinaga et al. [3] simulated the simplified sibilant
fricatives model to investigate the mechanisms of the production process in the vocal tract. Both the simplified
and realistic model are selected to be the validation by present framework. Fig. 13.23 is the instantaneous
velocity magnitude. Observably, the flow is accelerated at the narrow channel and the turbulence is caused
when the flow goes through the channel. As shown in Fig. 13.23, the sound source locates around the gap
before the exit, which is the high RMS value region. The fluctuation pressure out of the oral cavity is shownin Fig. 13.24, the frequency of this sound wave is about 5kHz. The SPL spectrum sampling at 10 cm far fromthe exit of oral cavity is shown as Fig. 13.25. The result simulated by current framework has good agreementwith the benchmark after 5kHz which is the first characteristic peak. Moreover, Fig. 13.26 is the instantaneous
velocity of realistic model. Clearly, the flow is speed-up between the articulators, and then products the soundby turbulence. Besides, the SPL spectrum shown in Fig. 13.27 has good tendency with the benchmark.
Present framework applied the implicit time scheme with immersed boundary method to achieves the highaccuracy and efficiency simulation in the computational acoustic of pronouncing sibilant fricatives. By validatingwith the benchmarks, the results have qualitatively good agreement. Moreover, the combination of immersed

144

CHAPTER 13. COMPLEX PHENOMENA UNIFIED SIMULATION RESEARCH TEAM

boundary method and high-performance computing reduces the computational resource drastically. It indicatesthat this framework is a promising tool for human phonation problem. About the future work, the investigationof how the shape effect the phonation will conduct for helping designing the denture.

Figure 13.22: Contours of fluc- Figure 13.23: (Top) instantaneous and (bottuation pressure for phase- and tom) RMS velocity magnitude.
spanwise-averaged flow fields.

Figure 13.24: Fluctuation pressure of simpli- Figure 13.25: SPL spectrum of simplified
fied model.
model.

Figure 13.26: Instantaneous velocity magniFigure 13.27:
tude of realistic model.
model.

13.3.3.4

SPL spectrum of realistic

References[1] Yokoyama, H., K. Kitamiya, and A. Iida. Physics of Fluids, 25.10, (2013), 106104.
[2] Nozaki K., Yoshinaga T., and Wada S., J Dent Res, 93(2), (2014), 207211.
[3] Yoshinaga, Tsukasa, Kazunori Nozaki, and Shigeo Wada. Physics of Fluids, 30.3, (2018), 035104.

145

13.3. RESEARCH RESULTS AND ACHIEVEMENTS

13.3.4

Speed-up and scale-up of fundamental methodology targeting post-K system

13.3.4.1

Performance optimization of incompressible solver of CUBE

This year, we evaluated and tuned the performance of the incompressible solver of CUBE on the K-computer.
The target problem is a simple cavity flow, which targets a standard numerical scheme used for empiricalanalysis (Multigrid + CG method). The serial tuning of the relevant part was performed by cutting out the
Poisson’s equation solver by the multigrid method, which is the most costly part in the incompressible solver(accounting for 91% of the total). Specifically, the following measures were applied.
1. The integer operation for address calculation has been reduced by making the array referenced
as a structure member directly referenced using a pointer.
2. Improved continuity of data access by swapping the dimensions of the array
3. The continuity of data access was improved by changing the array elements colored in red andblack by Red-Black ordering so that they are stored in separate arrays.
As a result, a performance improvement of 9.8 times was obtained. Next, the following two tunings wereapplied as parallel tuning.
1. In the packing process for neighbor communication, the cost imbalance between threads has
been improved by setting the division target of thread parallelization from each side of the
subdomain to the cube (collection of cells) included in each side.
2. By changing the OpenMP thread scheduling scheme from dynamic to static for the computational loop included in the CG method, the cost balance between threads is improved and theoverhead of dynamic scheduling is reduced.
As a result of these tunings (and the above-mentioned serial tuning), the performance of the entire program
was improved by 2.9 times at 8,192 nodes in the K computer. As a result, it can be said that the speedup
for creating results using ”Fugaku” has been completed. However, there are still unknown parts about theevaluation of scaling when the number of nodes is further increased and the expected increase in collective
communication cost, so we plan to deal with these in the result creation phase in the future.
13.3.4.2

Basic research for constructing a flow solution method using a neural network

In preparation for the coming era of computational science with AI initiatives, we are conducting basic researchon prediction technology based on neural networks, which is expected to be used in the future in this research
project. This year, we tried to reconstruct the flow field for the two-dimensional flow calculated by the unified
solution code CUBE based on the neural network (CNN) using the convolutional filter as a machine learningtechnology that can be used at present. The following network is constructed by following the method [1] of
replacing POD (Proper Orthogonal Decomposition) with a neural network for the flow field of the flow arounda cylinder.
!"#$% & ' (%&!#0$ + '

(%/ +
(/) +

!"#$) & ' ()&)*+ ' ,(- (/- +-.%Figure 13.28: Structure of neural network (partially modified figure in paper [1].
A neural network was trained for 10000 samples of snapshots as input, and the output was evaluated. As aresult, the mode decomposition equivalent to that of reference [1], that is, the flow field was decomposed into
two modes by a neural network, and the original flow field was reconstructed by synthesizing them, and the

146

CHAPTER 13. COMPLEX PHENOMENA UNIFIED SIMULATION RESEARCH TEAM!"#$%&'(")%#&%*+,-.(/#01234/2()%&"#$%&'(")
41'05%0(436"%0(2$#37

8'9(:6;(365()%"#/6"%<= (33#3

Figure 13.29: Example of flow field reconstruction results.

accuracy was equivalent to that of reference. It was confirmed that the flow field could be reconstructed. Figure13.29 shows an example of the results.
Currently, only 2D is handled due to computational resources and technical restrictions, but we plan toexpand this method in the future and evaluate the effectiveness of this method for 3D problems.

13.3.4.3

Unified solver extension using Adaptive Mesh Refinement (AMR)

Adaptive Mesh Refinement (AMR) has been found to be useful and effective for capturing fine structures influid analysis, such as where the gradient of a fluid state variable is steep, typiclly around a fluid-solid interface.
Since the calculation grid is dynamically subdivided only where necessary for the solution and shape of the flowthat changes quickly, it is not necessary to arrange a high-resolution calculation grid in the entire calculation
domain. Costs are significantly reduced, espacially in the unified multi-physics analysis such as the moving
boundary problem. However, the calculation load balancing for high parallel computing is well known as adifficult task to implement.
The CUBE developed in this research project already has a hierarchical data structure, and has a solution
for complex shapes, complex problems (heat, sound, chemical reaction, etc.) and moving boundaries. So, thebest approach to introduce the AMR method is to build an extension of the solver. To examine the directionof the approach, we conducted a feasibility study based on several well-known open source softwares, including
AMReX [2], which is being developed mainly by the Lawrence Berkeley National Laboratory in the UnitedStates. As a first step, we ran some basic simulations in search of new features suitable for our code. Figure13.30 shows an example of the calculation applied to the vortex movement problem. It can be seen that astime progresses, the subdivided regions of the grid also move. Early studies provided important insights ininvestigating the low implementation effort and missing features in the available code. In the future, we plan
to continue more detailed investigations and finally consider implementing it into CUBE.

Figure 13.30: An adaptive mesh refinement solution for a typical advection problem.

13.3. RESEARCH RESULTS AND ACHIEVEMENTS
13.3.4.4

147

References[1] T. Murata, K. Fukami, and K. Fukagata, “Nonlinear mode decomposition with convolutional neural networks
for fluid dynamics,” J. Fluid Mech. 882, A13 ,2020.
[2] AMReX-Codes, https://amrex-codes.github.io/, (cited 1st mar. 2020).

13.3.5Strategic development and deployment of industrial software for the construction of digital engineering test beds

13.3.5.1

Realistic automotive engine cycle simulation for the evaluation of the mixing efficiency

As part of the development of an automobile engine testbed, we performed a fuel-air mixing analysis thatassumed a practical usage scenario. The fuel-air mixing efficiency plays an important rule for automotiveengine design. Generally, through the observation of the fluid structures of the tumble and the swirl shown in
Figure 13.31, the efficiency can be evaluated.

Figure 13.31: The tumble and the swirl.
However, this kind of the turbulent structure is not easy to obtain from the experience. Therefore, the aimof this study is to simulate the cold flow inside the engine using CUBE and visualize the turbulent structure tocapture the tumble and swirl. Finally, the different shape of the piston is simulated to evaluate the performance
to show that this framework is suitable for the industrial applications. The governing equation and numericalmethod of this study is as follows.
∂U
∂F1
∂F2
∂F3
+
+
+
=0
∂t
∂x1
∂x2
∂x3
The quantities included in U and Fi areρ ρu1 
U = ρu2  ,
 ρu3 
ρe




ρui
 ρui u1 + P δi1 − µAi1





F =  ρui u2 + P δi2 − µAi2 ρui u3 + P δi3 − µAi3ρ(e + P )ui − µAij uj − k∂T /∂xi

where Aij = ∂ui /∂xj +∂uj /∂xi −2/3(5uδij ) and P is the pressure given by the ideal gas equation, P = ρRT .
The low-Mach-fix for Roe (LMRoe) scheme [1] with 5th order MSUCL reconstruction [2] is utilized in convection
term and second order central difference is used for the viscous term. For the moving geometry, the immersedboundary method (IBM) which has been implemented in CUBE [3,4] is adopted to the engine motion including
intake valve, exhaust valve and piston movement. First, the standard shape of the engine is simulated and theconservation is calculated to validate the practicality of the framework. The computational parameter and thephysical model are shown in Figure 13.32. Three different time steps are adopted.
The comparison of the results with the theoretical values from ideal gas state equation is shown in Fig.
13.33. According to the result, when the CF L < 10, the distribution of the density, pressure and temperature
has acceptable agreement with the theoretical result. The suggested reason for the discrepancy is that using
the IBM causes the leakage.
Next, the comparison of the different shape of the piston has been conducted. The computational parameterand the different shape of the piston are shown in Figure 13.35. The visualization result of the streamline withthe contour of the velocity magnitude is utilized to evaluate the magnitude of the tumble and swirl.

148

CHAPTER 13. COMPLEX PHENOMENA UNIFIED SIMULATION RESEARCH TEAM

Figure 13.32: The simulation conditions of the standard shape.

Figure 13.36 shows the results of the visualization result of the streamline with the contour of the velocity
magnitude. The result shows that with the convex at the left side of the piston in shape2, the velocity will bedecreased. On the other hand, shape1 and 3 still keep strong tumble. However, at the left side of the engine,
the tumble in shape3 becomes very weak. Overall speaking, shape1 obtain the strongest tumble.
Figure 13.37 shows the results at CA270 (Crank angle at 270 degrees). The result shows the same trend as
Fig. 5 that shape1 is the best shape. However, with the convex shape, the swirl will be enhanced because theflow field will not be symmetric anymore. The flow field in shape3 is not only in vertical direction but also inhorizontal direction so the swirl structure can be confirmed.
The results show that the current framework can be a useful tool for the engine design by giving a detailinsight of the flow mechanism. However, the calculation condition is still under the cold flow condition. In next
step, we are planning to model the combustion inside the engine and conduct direct evaluation of combustionenergy efficiency. And also we have started to provide this framework to the industrial user for the evaluationof the usability of the framework. We will continue to improve practicality through the consortium activities.

13.3. RESEARCH RESULTS AND ACHIEVEMENTS

Figure 13.33: The comparison of the density, pressure and temperature.

Figure 13.34: Engine cycle condition and piston shapes.

Figure 13.36: The streamline with the contour of the velocity magnitude at BDC.

149

150

CHAPTER 13. COMPLEX PHENOMENA UNIFIED SIMULATION RESEARCH TEAM

Figure 13.35: Engine simulation conditions.

Figure 13.37: The streamline with the contour of the velocity magnitude at CA270.
13.3.5.2

Numerical simulation of flow in a fuel-injector of an aircraft engine

As a first step of the development of the combustion system test bed, the flow simulation of the jet engine
combustor was performed. An investigation of grid sensitivity on local mesh refinement of numerical simulationof cold-flow in a fuel-injector of an aircraft engine was performed. In the case of an aircraft engine, a criticalvariable that governs emission, efficiency, and reliability is the equivalence ratio inside the combustion chamber.
This value is transported by swirling turbulent flow from fuel-injector to the combustion chamber. Therefore,prior to precise prediction and appropriate control of the combustion characteristics of the engine, as a first step,
the transport characteristics of flow through fuel-injector can be investigated through cold flow simulations. A
complete evaluation of flow characteristics of the engine requires numerical simulation with the entire annular
combustor geometry, which demands massive computational resources. Thus, an optimization strategy tominimize the computational cost through mesh optimization is warranted.
The cold flow simulations were carried out using CUBE, which is developed through this research. CUBE
is a multiphysics simulation framework designed for large scale industrial simulations. A hierarchical meshing151

13.3. RESEARCH RESULTS AND ACHIEVEMENTS

(a)

(b)

Figure 13.38: (a) The PIV measurement setup, (b) Numerical mesh and boundary condition

technique known as building cube method (BCM) is the foundation of CUBE over which numerical solvers
are built. Even with the use of BCM for meshing the fuel-injector geometry, using a same near wall spacing
mesh over the entire geometry results in excessively large numerical mesh. Consequently, relevant geometriccomponents of the injector that strongly influence the flow have to identified for mesh refinement while coarser
mesh can be used for rest of the geometry. Then it is important to investigate the dependence of flow on therefined mesh and to minimize the total mesh size while ensuring accuracy. The governing equations used in thisstudy are same as previous section. This is fully compressible flow approach. A 2D PIV measurement of air
flow in a fuel-injector was conducted by Kawasaki Heavy Industries. The detail flow field results were comparedwith the measurement result. The experimental setup and the measurement window are shown in Fig. 13.38.
In order to validate the accuracy, a very high-resolution mesh was used. The geometry used in the simulation
is highlighted as the fuel injector in Fig. 13.38. The boundary conditions and the mesh are also shown in
Fig. 13.38. The wall boundary condition, which means the boundary condition on the surface of geometry, wasset as adiabatic condition. The inflow mass flow rate was set as the same value of PIV measurement and theambient pressure 101,300 Pa, ambient temperature 288 K was set. The finest mesh resolution was 0.05 mm andthe number of calculation mesh was 500 million. The 2048 nodes of K-computer was used and the calculation
time was 48 hours to obtain the physical solution time of 2.7 s.

(a)

(b)

Figure 13.39: (a) Time averaged velocity magnitude high-resolution simulation, (b)Time averaged velocity
magnitude 2D PIV data
Velocity magnitude contours on a mid-sectional plane from the PIV measurement and the numerical simulation are shown in Fig. 13.39. Four distinct identifiable flow signatures, in the exit plane are identified in

152

CHAPTER 13. COMPLEX PHENOMENA UNIFIED SIMULATION RESEARCH TEAM

Fig. 13.40. The flow signatures are termed as A, B, C and D, in the PIV result, and A’, B’, C’ and D’ in thesimulation result. First, A (and A’) is the primary jet emanating out of the injector. A’ shows greater velocity
magnitude than A. In both cases, the flow gradually decelerates as it heads downstream. The flow signature B
is the formed by the flows from the pilot region and air-silt nozzle. It is observed that B is separated from A
and forms vortex. On the other hand, B’ does not form any vortex but just joins into A’, supplying momentum
to A’. It is hypothesized that intrinsically, B’ is just sucked into A’ because of the large velocity magnitudein the exit-vicinity region of A’. Meanwhile, C’ is more widely distributed although C is localized in narrow
region. It is regarded that the absence of vortex in B and the intensity difference primary jet (A-A’) may affectthe distribution. Furthermore, the velocity magnitude of the flow in C’ greater than the flow in C and it islikely that this is due to the difference in intensity of the primary jet, i.e. A-A’, which induces C-C’. Finally,
D’ stretches longer than D up to the exit of the injector. It can be considered that the greater intensity in theregion C’ and the absence of vortex in B’ assists the flow in D’ stretching it longer than its PIV counterpart D.

Figure 13.40: A comparision of the simulation and PIV results: 4 prominent features can be identified in eachcase.
To sum up, the overall difference of the velocity magnitude is due to that PIV data does not contain thedepth component. Meanwhile, it is regarded that the local differences are mainly affected by the difference
between B and B’ induced by overestimation of A’. Therefore, this is the sole qualitative difference between thesimulation and the PIV data. It can be inferred that the simulation qualitatively resolves the injector flow. In
the next step, we are planning to develop the fuel spraying model and the combustion model which is specificallydesigned for the real engine combustors. We will continue to develop the functionality to reproduce the realstate of engine combustion through the consortium activities.
13.3.5.3

Construction of surrogate model by aerodynamic multipurpose shape optimization
analysis and machine learning

In this study, we have developed a multipurpose optimization analysis framework using CUBE for the purpose of
efficiently optimizing aerodynamic performance at the initial design stages. This framework provides CHEETAH[5], a multipurpose evolution optimization software, in collaboration with the Japan Aerospace Exploration
Agency. CUBE has a data structure with a hierarchical Cartesian grid, and MeshWorks [6] is applied, which
is a shape morphing software by Detroit Engineered Products, Inc. It is a method that can be shared with asingle calculation grid for a plurality of changing vehicle shapes, and can avoid problems related to calculation
grid generation and morphing.
Continuing from last year, we performed genetic-evolutionary computation up to the 12th generation witha 4-objective function with the aim of obtaining an index of robustness against changes in the wind directionof sedan-type automobile models. The population size was 18, and a total of 252 cases of 36 + 18 + 18 × 12
generations including individuals in the L36 orthogonal table for generating the initial generation population.
And 2 wind directions for each needed to be calculated. Each case was calculated with 367 nodes of “Kcomputer” with 14 hours, and the total amount of calculation was about 2.6 million node-hours (including the13.3. RESEARCH RESULTS AND ACHIEVEMENTS

153

amount used last year). Figure 13.41 shows the Pareto solution for the difference between the aerodynamic
resistance value (Cd) and the resistance to the wind direction (∆Cd). The tendency of the objective function
value similar to the knowledge obtained in the actual design was reproduced, and the validity of the framework
could be confirmed.

Figure 13.41: The geometry deformation parameters (left) and optimized results of Pareto solution on 12th
generation (right).

In this physcal year, we tried to build a surrogate model by machine learning utilizing these results as
training-data. For neural networks constructed using Sony’s Neural Network Console [7], 18 individuals × 1
to 6 generations = 108 models were used for training-data, and 18 individuals × 7 to 12 generations = 108
models were used for evaluation-data. The input was 9 design variables, and Cd and ∆Cd were targets. Figure13.42 shows the aerodynamic drag obtained by calculation and the predicted result by machine learning. The
Cd value was predicted almost accurately, and it could be confirmed that the constructed neural network wasappropriate. It has also been confirmed that the Pareto solution based on this prediction result shows the sametendency as the solution obtained by the actual CFD calculation. This suggests that machine learning can bea powerful design tool in the future if appropriate design variables can be selected. However, there are somemodels whose prediction results are different with respect to the difference in aerodynamic drag, and there arealso subtle differences in the obtained sedan-type optimum body shape. We are proceeding with consideration
including investigation of the factors causing it. In the future, we plan to introduce design variables that
reproduce more detailed shapes and improve the degree of freedom, which will lead to the creation of results in“Fugaku”.

Figure 13.42: The prediction results using machine learning (left) drag (right) drag difference between 0 and -3
yawing angle.

154
13.3.5.4

CHAPTER 13. COMPLEX PHENOMENA UNIFIED SIMULATION RESEARCH TEAMWind resistant design of building structures – fluid structure coupled analysis of highrise towers

As an issue of the architectural consortium, verification analysis of a test bed based on a unique fluid structure
coupled analysis method was performed. In recent years, with the increase in the height of buildings and the
reduction in the weight of construction materials, there have been increasing cases of construction of high-risebuildings and large-scale structures having roofs using light materials. Since the design of these structures
requires low attenuation due to weight reduction, it is necessary to examine the wind resistance as well as theseismic resistance for earthquake. In particular, unstable aerodynamic vibrations caused by the interaction
between the vibrations of the building itself and the flow field around the building due to the turbulence of
the wind have a risk of causing a large amplitude leading to collapse of the building. Generally, a building isheavy and hard to be deformed as compared with a bridge or other structure. Therefore, conventionally, there
has been no problem with dynamic vibration, and it has been a problem to secure the building resistance to astatic wind load. However, with the emergence of such high-rise buildings, it has become important to study
the problem of vibration when the buildings are subject to strong winds. Thus, it is important to grasp the
structural characteristics in wind load evaluation. Figure 9 shows the mechanism of the generation of unstablevibration generated by the interaction between wind turbulence and the vibration of the building itself. Vigorous
research has been conducted on the unstable vibration phenomenon of buildings in the past, however there arestill many unknowns about the mechanism of occurrence and the influence of the interaction between them, andaerodynamic instability. It has been pointed out that when studying unstable phenomena such as vibration, thephenomena cannot be reproduced unless the structure and the surrounding flow field are treated as a strongly
coupled problem. There is few software that can achieve this. In this project, we have developed a unique
fluid-structure strong coupled unified solution that combines marker particles with the Euler-type structural
analysis method. In order to prove the validity of the method, detailed verification analysis was performed
using basic tower structure.

Figure 13.43: Unstable vibration phenomenon of high-rise building.
In order to verify the validity of this method, first, a fluid-structure interaction analysis with a rigid structure
was performed. The structural shape is a quadrangular prism of 0.1 × 0.1 × 0.1m. We compare it with theresults of previous experiments. The calculation conditions, analysis mesh and inflow fluctuation wind are as
follows.
Figure 13.48 shows a comparison with the results of previous wind tunnel experiments. The black line
represents the result of the calculation by this method, and the blue line represents the result of the previous
experiment, and it can be confirmed that the values agreed well.
Considering the introduction of this method to architectural structures, numerical stability was examinedby performing analysis with various conditions of prismatic structures. Four cases were prepared in which the
Reynolds number was fixed at 10, the mass density difference between the solid and the fluid, and the shear
modulus of the solid has been changed. The analysis results are shown in Fig. 13.49-13.50. From these results,
it was found that numerical instability occurred near the prisms in all models due to the numerical calculation
method of solid strain. In the next step, we plan to introduce a method that can stably calculate the solid strain
without depending on the velocity field. We will continue to develop the functionality to produce the practicallyuseful framework for the strong-coupling fluid-structure interaction through the consortium activities.

13.3. RESEARCH RESULTS AND ACHIEVEMENTS

155

Figure 13.44: Calculation conditions.

Figure 13.45: CUBE mesh.

Figure 13.50: Comparison of calculation staFigure 13.49: Mises stress distribution of bility of each case.
CASE-D.

Finally in this fiscal year, we conducted a simulation of wind resistance analysis of the dome roof of a stadium
structure as an application to practical problems. The flow condition is 10 m/s of air wind, the Reynolds number
is approximately 273 million, and the soft material property is set as 33,300 Pa of elastic modulus, 300kg/m3
of the mass density of solid. The number of the cell was approximately 17 million, and the calculation time was
10 hours using 256 nodes of K-computer. Figure 13.51 shows a typical flow field coupled with the soft material
deformation of the stadium roof. We cannot say anything about the validity of this result or its physics yet.
However, this result may be the world’s first example of solving FSI problems at a high Reynolds number overhundred million by an Eulerian FSI method. At least, we did not observe the unphysical deformation of theroof region (sponge) induced by wind force. In the next step, we plan to increase the resolution of calculation
mesh to try to capture smaller flow structures and deformation modes. We will continue to contribute to help
users to apply this method into the practical cases through the consortium activities.

156

CHAPTER 13. COMPLEX PHENOMENA UNIFIED SIMULATION RESEARCH TEAM

Figure 13.46: Inflow fluctuating wind profile.

Figure 13.47: X-direction velocity distribution in vertical section of prism.

Figure 13.48: Horizontal mean wind pressure coefficient distribution of prism: comparison with previous experimental results.

Figure 13.51: Example of flow and stress results on deformable roof analysis for a stadium structure.
13.3.5.5

Topology optimization using Eulerian structure analysis

Topology optimization [8] is a method to determine the arrangement of materials and components, including
the voids, to maximize or minimize the objective function using a limited material volume in the design space.
In recent years, industrial applications of topology optimization have been increasing as it provides effective
suggestions to designers. On the other hand, when considering the application to a real product, the volume of
the actual product material is usually less than 1% of the design area expected by the designer. In addition tothis, the preparation of design-domain shape modeling, which includes many cases and trials in the upstreamof the design process, requires a lot of time and effort. In this study, we are developing an original unified fluidstructure analysis method based on Eulerian methods combined with the Lagrangian marker particle method
to solve a large-scale structural analysis problem which is difficult to be solved by conventional methods. The

157

13.3. RESEARCH RESULTS AND ACHIEVEMENTS

marker particle method is one of the most promising methods to solve the topology optimization challenges
mentioned above, because it is easy to achieve high parallelization efficiency and to reduce the time and effortof mesh generation in a massively parallel computing environment. In this year, we have developed a topology
optimization framework using Eulerian structural analysis method. To verify the validity of the topology
optimization, the optimization of the maximum amount of energy absorbed by a two-dimensional structurewith a circular hole in the center of the design region under enforced velocity loading in the direction of shear
was carried out. The results are comparable to those obtained by the finite element method, even when CUBE
is used. In addition, the quantitative value of the objective function also approaches that of FEM by increasing
the spatial resolution.

(a)

(b)

Figure 13.52: Optimal shape as a result of topology optimization to maximize absorption energy by shear
enforced velocity load (left), and its objective function history (right).

Figure 13.53: Design area which surrounds a conventional white body frame (left) and optimal shape as a resultof topology optimization subject to torsional enforced velocity load (right).
Using this framework, we performed optimization trials for a car body frame with 10 million design variables,
which is somewhat difficult to do with existing topology software (Fig. 2). A design domain was set up over thewhite body frame of an existing published body model (TOYOTA, Yaris 2010, [9]), and forced-velocity loads
were applied to the design domain to maximize the amount of energy absorbed by the body torsion. As a result,a frame structure close to the intended white body frame was obtained. In the future, we will continue to carryout various basic verification examples and optimization verifications for a large number of load cases that are

158

CHAPTER 13. COMPLEX PHENOMENA UNIFIED SIMULATION RESEARCH TEAM

assumed to be real problems, as well as to improve the optimization performance through collaboration withindustry.
13.3.5.6

Implementation of aerodynamic body motion coupled analysis framework

In recent years, it has become clear that the evaluation of unsteady aerodynamic force that reproduced in theactual running state is indispensable for the aerodynamics parts design, since it affects the evaluation of steeringstability of automobiles. In general, it is difficult to accurately capture the flow field around a vehicle during
actual driving motion and to measure changes in vehicle attitude. And the actual driving test is affected byatmospheric disturbances, so the same situation should be reproduced. Additionaly, it is difficult to evaluate thesteering stability by quantitative evaluation depending on the driver’s experience. Therefore, in this research,
we constructed an analysis framework that can reproduce the actual driving state by combining aerodynamic
analysis with vehicle motion on a supercomputer, which includes suspension, steering mechanism and driver’s
operation reaction, etc.
This framework adopts the implementation of the open source multibody dynamics analysis software, named
Project Chrono [10], into the CUBE. It is possible to use commercial vehicle motion analysis software, but inorder to operate it on a supercomputer, we have experienced problems such as binary support, license handling,and external communication of calculation nodes. So an open source framework is powerful solution, thesignificance of being able to do it is high rather than other researchers think. In this year, we succeeded inconstructing this framework and implemented a practical example of vehicle motion analysis.
A hatchback type real vehicle shape model was used for the analysis. The model also reproduces detailed
shapes such as underfloor and auxiliary equipment inside the engine room. We are comparing steering stability
with and without aero parts mounted on the side of the vehicle. In this case, the suspension and steering
mechanism of the same type as the actual vehicle could be used. Assuming a lane change movement to theovertaking lane-chage maneuvering on the highway, the steering stability was evaluated for the movement to
the right side 3 m at 100 km/h.
Figure 13.54 shows an example of visualization results of flow field and motion analysis. A correlation
has been shown between the sensory evaluation of steering stability and the roll motion of the vehicle. Whenevaluated by focusing on the roll motion and the steering input by the driver (Fig. 13.55), both the angularvelocity and the angular jerk (time change of acceleration) showed a decrease in peak and amplitude when theroll motion in the first half converged. Results consistent with the report were obtained in the driver’s sensory
test. In the future, we plan to make improvements to create results on “Fugaku” while conducting higher
resolution analysis and improving versatility.

Figure 13.54: Analysis example of lane change maneuver (upper: CUBE result, lower: motion analysis visualization result).

13.4. SCHEDULE AND FUTURE PLAN

159

Figure 13.55: Results of roll angular velocity of vehicle motion.
13.3.5.7

References[1] F. Rieper, “A low-Mach number fix for Roe’s approximate Riemann solver,” J. Comput. Phys., 230, pp.
52635287, (2011).
[2] K.H. Kim, C. Kim, “Accurate, efficient and monotonic numerical methods for multi-dimensional compressible
flows Part II: Multi-dimensional limiting process,” J. Comput. Phys., 208, pp.570615, (2005).
[3] C.-G. Li, M. Tsubokura, R. Bale, “Framework for simulation of natural convection in practical applications,”
Int. Commun. Heat Mass Transf., 75, pp.5258, (2016).
[4] N. Jansson, R. Bale, K. Onishi and M. Tsubokura, CUBE, “A scalable framework for large-scale industrial
simulations, International,” Journal of High Performance Computing Applications, Vol.33, pp.678-698, (2018).
[5] A. L. Jaimes, A. Oyama, and K. Fujii, “A ranking method based on two preference criteria: Chebyshev
function and e-indicator,” Proceedings of 2015 IEEE Congress on Evolutionary Computation, (2015).
[6] MeshWorks WEB site, https://www.depusa.com/meshworks/, [cited 1 Mar 2020].
[7] Neural Network Console WEB site, https://dl.sony.com/ja/, [cited 1 Mar 2020].
[8] M.P. Bendsøe, N. Kikuchi: “Generating optimal topologies in structural design using a homogenization
method,” Comput.Methods in Appl.Mech.Eng. Vol 71, No. 2, pp. 197 224 (1998).
[9] Center for Collision Safety and Analysis: “2010 Toyota Yaris Coarse Finite Element Model Version 1l”, 2016
Dec. URL: https://www.ccsa.gmu.edu/models/2010-toyota-yaris/[10] Project Chrono WEB site, https://projectchrono.org/, [cited 1 Mar 2020].

13.4

Schedule and Future Plan

In the research and development of unified simulation framework ’CUBE’ from 2012, the following goals has
been achieved:
1. Construction and development of the simulation technology for bringing out the performance of Kcomputer
2. Validation and practical usage of industrial applications such as vehicle aerodynamics, aeroacoustics,
unsteady motion aerodynamics, engine combustion, city area wind environment analysys, sport CFD, andstructure analysis
3. Preparation of the simulation technologies of HPC toward EXA-scale
In the long-term objectives, the following target is considered:
1. Establishment of the research and development center for industrial simulation technology
2. Contribution to computer science by expanding the developed simulation technology to different fields
In terms of using Cartesian grids, the problems that the analysis accuracy in the vicinity of the wall surface
decreases due to the limitation of immersed boundary method is becoming clear. In the future, through the

160

CHAPTER 13. COMPLEX PHENOMENA UNIFIED SIMULATION RESEARCH TEAM

deployment of CUBE to the industrial applications, the introduction of higher order schemes, wall modelingand advance sophisticated research of immersed boundary method are planning.
In terms of structure analysis, to avoid this numerical oscillation even when the spatial resolution is low, weare going to introduce a calculation method for solid deformation tensors that does not require velocity gradientcalculation and we continue to work on development to improve analysis accuracy.
In terms of using Hierarchical grid system, the problems that the interface communication in each halo of
cubes will spoil the performance on future supercomputers is becoming clear. In the future, the enhancement
to the current communication / calculation overwrapping technique, or tuning for the interface communication,
introduction of adaptive refinement scheme are planning.
From collaboratetive companies, including companies through the activities of the “Consortium for Next
Generation Automotive CAE using HPC” and “Consortium for Next Generation Combustion CAE using HPC”
organized by RIKEN, voices of surprises and expectation for software development have been received. We areplanning to continue development for practical application in near future.

13.5

Publications

13.5.1

Articles/Journal[1] Rahul Bale, Neelesh A. Patankar, Niclas Jansson, Keiji Onishi, and Makoto Tsubokura, “Stencil Penalty
approach based constraint immersed boundary method,” Computers & Fluids, Vol.200(30), 104457, (2020).
[2] Koji Nishiguchi, Tokimasa Shimada, Masafumi Otaka, Shigenobu Okazawa, Makoto Tsubokura, “Eulerianfinite volume formulation using Lagrangian marker particles for compressible solid analysis,” Journal of Japan
Society of Civil Engineers, Ser. A2 (Applied Mechanics (AM)), 75(2), I 237-I 248, (2019).
[3] Keiji Onishi, Makoto Tsubokura, “Topology-free immersed boundary method for incompressible turbulence
flows: An aerodynamic simulation for ’dirty’ CAD geometry”, arXiv preprint, arXiv:2002.06206, (2020)

13.5.2

Conference Papers[4] Keiji Onishi, Rahul Bale, Makoto Tsubokura, “Assessment of Rotating Wheel Vehicle Aerodynamics Simulation using Cartesian Grid Method and Open-grill Full Vehicle Models,” SAE Technical Paper 2019-01-0660,(2019).
[5] Onishi Keiji, Jansson Niclas, Tsubokura Makoto, “In-situ visualization on large scale vehicle aerodynamics
simulation,” In: Proceedings of the 24th Conference on Computational Engineering and Science, A-05, (2019).
[6] Rahul Bale, Kazuto Ando, Yuki Hayashi, Keiji Onishi, and Makoto Tsubokura, “Analysis of underwater
undulatory swimming using self-propulsion simulations,” In: Proceedings of the 24th Conference on Computational Engineering and Science, (2019).
[7] Koji Nishiguchi, Tokimasa Shimada, Masafumi Otaka, Shigenobu Okazawa, Makoto Tsubokura, “Eulerianfinite volume formulation using particle-in-cell method for large-scale parallel simulation of complex structures,”
Proceedings of ECCOMAS MSF 2019 THEMATIC CONFERENCE, 4 pages, Sarajevo(Bosnia-Herzegovina),(2019).
[8] Koji Nishiguchi, Rahul Bale, Tokimasa Shimada, Masafumi Otaka, Hirofumi Sugiyama, Shigenobu Okazawa,
Makoto Tsubokura, “Eulerian finite volume formulation using Lagrangian maker particles for compressible solid
analysis,” In: Proceedings of the 22nd JSCE Applied Mechanics Symposium, C000134, (2019).
[9] Koji Nishiguchi, Tokimasa Shimada, Masafumi Otaka, Shigenobu Okazawa, Makoto Tsubokura, “Structure
analysis with particle-in-cell method using hierarchical Cartesian mesh,” In: Proceedings of the Conference on
Computational Engineering and Science 24, A-05-02, (2019).
[10] Tokimasa Shimada, Koji Nishiguchi, Shigenobu Okazawa, Makoto Tsubokura, “Eulerian Unified Formulation Using Lagrangian Marker Particles for Fluid-Structure Interaction Problems,” In: Proceedings of the 33th
Computational Fluid Dynamics Symposium, E11-4, (2019).
[11] Hiromasa Kajimoto, Shunpei Koyama, Rahul Bale, Keizo Yamamoto, and Makoto Tsubokura, “Aerodynamics simulation framework for Ski-jumping take-off and its application to various jumpers,” In: Proceedingsof the 33th Computational Fluid Dynamics Symposium, D01-4, (2019).
[12] Yotaro Nomoto, Rahul Bale, Keiji Onishi and Makoto Tsubokura, “Construction of a Car AerodynamicMotion Coupled Analysis Framework Using Hierarchical Orthogonal Grid,” In: Proceedings of the 33th Computational Fluid Dynamics Symposium, B07-1, (2019).
[13] Younghwa Cho, Rahul Bale, Chung-Gang Li, Tomonari Satoh, Ryusuke Matsuyama, Takeo Oda, Makoto

13.5. PUBLICATIONS

161

Tsubokura, Nobuyuki Oshima, “Numerical simulation of flow in a fuel-injector of an aircraft engine,” In: Proceedings of the 33th Computational Fluid Dynamics Symposium, (2019).
[14] Tetsuro TAMURA, Yong CAO, Yuki NAGAO, Hidenori KAWAI, Koji NISHIGUCHI, Makoto TSUBOKURA, “High performance computing of wind-induced pressures on various domes under real conditions,”
Proceedings of the IASS Annual Symposium 2019 (FORM and FORCE 2019, Barcelona, Spain), pp.1185-1192,(2019).

13.5.3

Posters[15] Kazuto Ando, Rahul Bale, Keiji Onishi, Kiyoshi Kumahata, Kazuo Minami, Makoto Tsubokura, “Optimizing Multigrid Poisson Solver of Cartesian CFD code,” SC19 - International Conference for High Performance
Computing, Networking, Storage and Analysis, Denver(CO), (2019).
[16] Onishi Keiji, Bale Rahule, Tsubokura Makoto, “Research and development of automotive aerodynamic
simulation method for WLTP,” Japan Society of Mechanical Engineers, the 97th Fluid Engineering Division
Conference, (2019).
[17] Bale Rahul, Wang Wei-Hsiang, Li Chung Gang, Onishi Keiji, Uchida Kenji, Hidefumi Fujimoto, Tsubokura
Makoto, “Numerical Simulation of Combustion in a Rapid Compression Machine (RCM),” Japan Society of
Mechanical Engineers, the 97th Fluid Engineering Division Conference, (2019).
[18] Ando Yasunori, Li Chung Gang, Onishi Keiji, Tsubokura Makoto, “The Prediction of Narrowband Noise
Emitted from Real Vehicle Shape,” Japan Society of Mechanical Engineers, the 97th Fluid Engineering Division
Conference, (2019).

13.5.4Invited Talks

13.5.5

Oral Talks[19] Keiji Onishi, Yasunori Ando, Kosuke Nakasato, Makoto Tsubokura, “Assessment of Rotating Wheel Vehicle
Aerodynamics Simulation using Cartesian Grid Method and Open-grill Full Vehicle Models,” WCX19: SAE
World Congress Experience, Detroit (MI), (2019).
[20] Keiji Onishi, Makoto Tsubokura, “Toward Realization of Real-Time Vehicle Aerodynamics Simulation,”
ASME-JSME-KSME Joint Fluids Engineering Conference (AJK Fluids2019), San Francisco, CA, (2019).
[21] Koji Nishiguchi, Rahul Bale, Tokimasa Shimada, Masafumi Otaka, Shigenobu Okazawa and Makoto Tsubokura, “Eulerian Finite Volume Formulation using Lagrangian Marker Particles for Deformable Solid-Fluid
Interaction Problems,” Proceedings of the 7th Asia-Pacific Congress on Computational Mechanics (APCOM,
Taipei), (2019).
[22] Koji Nishiguchi, Tokimasa Shimada, Rahul Bale, Shigenobu Okazawa, Makoto Tsubokura, “Particle-In-Cell
Method using Hierarchical Cartesian Mesh for Deformable Solid-Fluid Interaction Problems,” Proceedings of
VIII International Conference on Coupled Problems in Science and Engineering, 1 page, Sitges (Spain), (2019).
[23] Shunpei Oyama, Keizo Yamamoto, Hiromasa Kajimoto, Rahul Bale, Jun Ikeda, Makoto Tsubokura, ”Unsteady aerodynamics simulation for the influence of difference in posture change on aerodynamics and flight
distance in ski jump take-off motion,” JSME Symposium : Sports engineering and Human Dynamics, B-2,
(2019).
[24] HsuehJui Lu, ChungGang Li, Akiyoshi Iida, Tsukasa Yoshinaga, Kazunori Nozaki, Makoto Tsubokura,“A framework for simulation of sibilant fricatives using implicit compressible flow solver,” The 72nd Annual
Meeting of the American Physical Society’s Division of Fluid Dynamics (APSDFD)(Nov 23-26 2019, Seatle),
C29.00005, (2019).
[25] HsuehJui Lu, ChungGang Li, Makoto Tsubokura, “The Simulation of Flow and Acoustics for Human Phonation System Using Implicit Compressible Flow Solver,” International Conference on Flow Dynamics(ICFD)(Nov
6-8 2019, Sendai), OS12-8, (2019).

13.5.6

Software13.5.7

Patents