Vol.4 ,No.1, May 2017

Transactions on

GIGAKUSpecial issue of
The 5th International GIGAKU
Conference in Nagaoka(IGCN), 2016

GIGAKU Press
Nagaoka University of Technology

Transactions on GIGAKU: Scope and Policy
Nagaoka University of Technology publishes an online, open access journal titled “Transactions on
GIGAKU”, which is focused on the science and technology related to GIGAKU*. The mission of thisjournal is to spread out the concept of GIGAKU and the fruits of GIGAKU to the global world and to
be a strong network for innovations in science and technology and for development of next
generations of high-level human resources. This journal, therefore, covers research and education
activities related to GIGAKU in broad areas.
*

See ‘What is GIGAKU?’ below.

‘What is GIGAKU?’GIGAKU is a term composed of two Japanese word-roots; GI and GAKU. The word GI〔技〕
literally stands for all kinds of arts and technology, and GAKU 〔学〕 stands for scientificdisciplines in general when used as a suffix.
The term was originally coined to describe the fundamental philosophy of education and research of
Nagaoka University of Technology (NUT) when it was established in 1976. Through this term thefounders of NUT intended to express their recognition that all technical challenges in the real world
require a scientific approach. And NUT has a relentlessly pursued GIGAKU since then.
Thirty-five years have passed and all surrounding conditions have changed dramatically during those
years. We are witnessing rapidly globalizing economics and huge scale changes in demographic,
industrial and employment structures. All those changes seem to necessitate the further evolution of
GIGAKU. In response to this, NUT recently announced its new “Growth Plan” and a reneweddefinition of the term is given;
GIGAKU is a science of technologies, which gives us an angle to analyze and reinterpret diverse
technical processes and objects and thus helps us to advance technologies forward. By employing abroad range of knowledge about science and engineering, management, safety, information
technology and life sciences, GIGAKU provides us with workable solution and induces future
innovations.

July 2012
Koichi Niihara, President of Nagaoka University of Technology

1Editors (2017 to 2018):
Editor-in-chief:
Prof. Yasuhiro Wada
Head of Library, Nagaoka University of Technology, Nagaoka Japan
Co-Editors:
Prof. Nobuhiko Azuma
President of Nagaoka University of Technology, Nagaoka Japan
Prof. Yoshiki Mikami
Nagaoka University of Technology, Nagaoka Japan
Prof. Yasuo Utsumi
Sendai National College of Technology, Sendai Japan
Prof. K. Bhanu Sankara Rao
Mahatma Gandhi Institute of Technology, Gandipet, Hyderabad India
Prof. Paritud Bhandhubanyong
Panyapiwat Institute of Management (PIM), Nonthaburi, Thailand
Prof. Takumi Fujiwara
Tohoku University, Sendai, Japan
Submission and Publishing“Transactions on GIGAKU” publishes original papers focusing on education, science,and technology related to GIGAKU. Proposals for special issues aligned withconferences are strongly encouraged*, and guest editors are welcomed. Papers
submitted are peer-reviewed before publication. Publication charge is free for all
submitted papers.
* Please contact to Editor-in chief, ywada@vos.nagaokaut.ac.jp
Copyright
The copyright of the papers published in Transactions on GIGAKU is transferred to
Nagaoka University of Technology. (See our copyright transfer form for details.)
Contactjoho-kiban@jcom.nagaokaut.ac.jp(Revised May 2017)2

Editorial board of special issue of The 5th International GIGAKU
Conference in Nagaoka (IGCN), Vol.4, No.1. of Transactions on
GIGAKU
EditorDr. Masato Aketagawa, Professor, Department of Mechanical,, Nagaoka University ofTechnology, Japan

Co-Editors
Dr. Shigeru Nagasawa Professor, Department of Mechanical, Nagaoka University ofTechnology, Japan
Dr. Takashi Yukawa, Professor, Department of Management and Information Systems
Engineering, Nagaoka University of Technology, Japan
Dr. Tomoichiro Okamoto, Associate Professor, Department of Electrical, Electronics andInformation Engineering, Nagaoka University of Technology, Japan
Dr. Yukio Miyashita Associate Professor, Department of Mechanical, Nagaoka
University of Technology, Japan
Dr. Yuichi Otsuka, Associate Professor, Department of System Safety, Nagaoka
University of Technology, Japan3

Transactions on GIGAKU
Volume 4, No. 1, May 2017
The 5th International GIGAKU Conference in Nagaoka (IGCN)(October 6-7, 2016, Nagaoka University of Technology, Nagaoka, Japan)
The IGCN is designed and organized to provide a cross-border, cross-sector,
cross-disciplinary forum for those researchers, educators, and industrial leaders who are
creating and practicing GIGAKU in various technology domains in various countries.
Some papers presented in the IGCN are published in this first volume of Transactions
on GIGAKU.
Table of contents
1. Pages 04001/1-8
Shape identification of cavity in structures using thermal testing data based on adjoint
variable and finite element methods
Kotaro Maruoka, Takahiko Kurahashi, Tetsuro Iyama
2. Pages 04002/1-7
Optimizing boundary surface of ceramics die in deep drawing process
Hanh C. Nguyen, Shigeru Nagasawa, Kensei Kaneko
3. Pages 04003/1-9
Bending behaviour of polypropylene sheet subjected to two-line wedge indentation
Akihiro Yamamoto, Satawat Singprayoon, Shigeru Nagasawa, Seksan Chaijit
4. Pages 04004/1-6
Improvement to Reduce Delay Time in Transportation to Disaster Area using Similarity
to Electric Circuit
Chata Yorvarak, Thanat Rungsirathana, Kazumasa Takahashi,
Toru Sasaki, Takashi Kikuchi, and Nob. Harada
5. Pages 04005/1-7
An Observation on Potential Alkali Silica Reactivity of Natural Aggregate in Myanmar
Pyae Phyo Kyaw, Thynn Thynn Htut
6. Pages 04006/1-8
Parallel implementation of Entropic lattice Boltzmann method for flow past a circularcylinder at high Reynolds number
Ayurzana Badarch, KhenmedekhLochin, Hosoyamada Tokuzo
7. Pages 04007/1-7
A Study on the Relationship between the Mindset of Taxi Drivers and the Attribution of
Accident Causes and Driving Behavior
Makiko Okamoto, Katsuko Nakahira,T4

8. Pages 04008/1-6
Engineering Design Exercise (EDE) for Faculty Development in National Institute ofTechnology, Nagaoka College
Rumi Shindo Togashi , Taku Kiryu , Shigehiro Toyama , Yoshinori Tokoi ,Tetsuro Iyama , Fujio Ikeda , Yuki Murakami , Shin-ichi Akazawa ,Yasuko Tsuchida5Transactions on GIGAKU 4(1)(2017)04001/1-8

Shape identification of cavity in structures using thermal testing data
based on adjoint variable and finite element methods
Kotaro Maruoka1,*, Takahiko Kurahashi2**, Tetsuro Iyama3
Nagaoka University of Technology, 1603-1 Kamitomioka-cho, Nagaoka 940-2188, Japan,
2)
Department of Mechanical Engineering, Nagaoka University of Technology,1603-1 Kamitomioka-cho, Nagaoka 940-2188, Japan3)
Department of Mechanical Engineering, National Institute of Technology, Nagaoka college,
888 Nishikatakai, Nagaoka 940-8532, Japan1)*E-mail: s133075@stn.nagaokaut.ac.jp**E-mail:kurahashi@mech.nagaokaut.ac.jp

In this study, we introduced the smoothing techniques into thethree-dimensional shape identification analysis of the cavity shape in test
piece based on thermal data and compared diference between computational
results employing each techniques. As a result, it is found that appropriatesmoothing technique on the shape identification analysis varies according toconditions in target shape. Gaussian filter is appropriate in a case where theshape of cavity is identified as including smooother surface from the aspectof performance function. On the other hand, in the case that the cavity surfaceof the target shape has a dent, traction method has efficacy for bringing thecomputational model close to target shape.

1. Introduction
In the field of the shape optimization, there is the problem that the calculated sensitivities for theobjective function become oscillation. Due to this problem, optimum shape is often obtained withconcavo-convex surface. This problem is also seen in the shape identification analysis proposed by
Kurahashi and Oshita[1], based on the adjoint variable method and the finite element method. In thisanalysis, identification of cavity shape is carried out by using temperature history on the thermal
testing method. In previous study, Kurahashi et al. applied shape identification analysis to the cavity
shape in test piece created using a 3D printer and introduced Gaussian filter, a kind of smoothing
process, to modify the oscillation of the gradient vector[2]. As a result, identified shape close to target
shape was obtained. However, appropriate condition for smoothing on the shape identification has
not been investigated. In this study, we carry out shape identification analysis with two smoothing
techniques, Gaussian filter and traction method[2]. The computational results are compared tocharacterize each processes on some computational conditions.

2. Shape identification analysis
The computation approach in shape identification analysis is shown in this section by reference
to previous research[3]. State equation is given as heat transfer equation (See Eq. (1)).
ߩܿ߶̇ − ߢ߶,௜௜ = 0(1)

For the heat transfer equation, the initial condition and boundary condition are defined as in Eq. (2).
෠⎧ ߶ = ߶(‫ݐ‬଴ ) in Ω
⎪ ߶ = ߶෠ on Γଵ⎨−ߢ߶,௜ ݊,௜ = ℎଵ ൫߶ − ߶௜௡௙ ൯ on⎪⎩−ߢ߶,௜ ݊,௜ = ℎଶ (߶ − ߶௖௔௩ ) on

1

Γଶ
Γଷ

(2)

Transactions on GIGAKU 4(1)(2017)04001/1-8

where ρ, c, κ, ni , h1, h2, ϕinf and ϕcav denote density, specific heat, thermal conductivity, unit normal
vector, the convective heat transfer coefficient of the space inside the cavity, the convective heattransfer coefficient of the outside surface, ambient temperature and the temperature inside the cavity.
Γ1 denotes the lower surface, Γ2 denotes the outside surfaces, and Γ3 denotes the surface of the cavity(See Fig. 1).

Fig. 1 Boundary definition of the test piece
The performance function is given as Eq. (3) to evaluate residual error between observational dataand computational temperature1 ௧೑‫ = ܬ‬න න ܴ(߶ − ߶௢௕௦ )ଶ ݀Ω ݀‫ݐ‬
2 ௧బ ஐ

(3)

From the state equation and the performance function, following Lagrangian function is obtained.
௧೑‫ ܬ = ∗ܬ‬+ න න ߣ(ߩܿ߶̇ − ߢ߶,௜௜ )݀Ω ݀‫ݐ‬
௧బ

(4)

ஐ

This function indicates performance function added constrained conditions, i.e., the state equation,
the boundary conditions and the initial condition. The first variation, which is a necessary condition
for minimizing the Lagrangian function, is shown in Eq. (5).
ߜ‫= ∗ܬ‬

߲‫∗ܬ‬
߲‫∗ܬ‬
߲‫∗ܬ‬
߲‫∗ܬ‬
ߜߣ +
ߜ߶̇ +
ߜ߶ +
ߜ‫ݔ‬
̇
߲ߣ
߲߶
߲‫ݔ‬௜ ௜߲߶

(5)where xi indicates coordinates on the design variable. The stationary condition is given as a condition
under which each term equals zero. Then final term in Eq. (5) indicates the gradient vector withrespect to the coordinates on the cavity surface. The cavity shape is updated with the gradientvector, and followed by recalculation of the state equation. The procedure for inverse analysis is
described below.
(1) Set of initial coordinates xi(0) and convergence criterion ε.
(2) Computation of state equation.
(3) Computation of performance function J(0).
(4) Computation of adjoint equation.
డ௃∗(೗)(5) Update of coordinate xi(l) based on the computed gradient vectorቄ డ௫ ቅ.
೔(l is the number of iterations, W is weighting parameter.)
(௟ାଵ)

ቄ‫ݔ‬௜

()

(௟)

ቅ = ቄ‫ݔ‬௜ ቅ − [ܹ (௟) ]ିଵ ቊ()߲‫(∗ܬ‬௟)ቋ߲‫ݔ‬௜(6) Checking for convergence; if (‫ ݈ ܬ‬− ‫݈ ܬ‬+1 )/‫(ܬ‬0) < ε, then stop, else go to step 7.
(7) Computation of state equation.
(8) Computation of performance function J(l+1)
2Transactions on GIGAKU 4(1)(2017)04001/1-8(9) Updating of weighting parameter W(l)I. if J(l+1) < J(l), then W(l) = W(l) × 0.9, and go to step 4
II. if J(l+1) > J(l), then W(l) = W(l) × 2.0, and go to step 5
3. Computational conditions
In previous study, thermal testing was carried out by using real test piece made of ABS. Fig. 2
shows detail of the test piece, and Fig. 3 shows finite element model. Considering the fact that it iseasy to make the test piece by the 3D printer, the size and the shape of cavity is decided. The
observation point is positioned on the center of the top surface. We describe the computational
condition that the cavity shape shown in Fig.2 is set to target shape ”Case A”. Thermal conditions forthe test piece and computational conditions are referred to previous study and given as shown in
Tabs. 1 and 2. In addition, using another test piece incorporating the cavity with a dent, the heattransfer analysis is carried out, and the temperature at the observation point is obtained. The drawingof the model and finite element model are shown in Figs. 4 and 5. The computational condition
setting target shape as the cavity shape shown in Fig.4 is described “Case B”. In computation ofshape identification, the initial shape model is given as incorporating a cavity thickness of 13 mm.

Density ρ

Table 1 Thermal properties of test pieceSpecific heat c Thermal conductivity κ
Heat transfer coefficient
Outside surface

3

ABS

2

Cavity[g/cm ]

[J/kgK][W/mK][W/m K][W/m2K]

0.75

1386

0.135410.0

7.0

Table 2 Computational conditions
Total number of nodes
1460
Total number of elements

5552

Total time tmax, sec

600

Time increment Δt, sec

10Time steps

60

Convergence criterion ε10-6Fig. 2 Drawing of the test piece (1)Fig. 3 Finite element model of the target shape (1)3Transactions on GIGAKU 4(1)(2017)04001/1-8

Fig. 4 Drawing of the test piece (2)Fig. 5 Finite element model of the target shape (2)4. Shape identification analysis introduced smoothing method
Shape identification analysis under the computational condition ”Case A” is carried out. First,
numerical computation without the smoothing is carried out. The identified shape of the cavity isshown in Fig. 6. The cavity surface is identified convexly. In addition, normalized performance
function decreases to about 0.55 (See Fig. 7). Thus the identified shape similar to target shape cannot
be obtained. Next we applied Gaussian filter to gradient vector. Gaussian filter is one of the digital
filter used in the field of graphics[4]. In this smoothing method, smoothed values are determined of
all the nodes on target domain and weight parameter shown in Eq. (6)‫ݔ(ܩ‬, ‫= )ݕ‬ି௫ మ ା௬ మ
1
ି
ଶఙ మ
݁
2ߨߪ ଶ

(6)

where x and y indicate the x-axial distance and the y-axial distance between the target node andreferenced node. In this study, the weighting parameter for node m to node n is defined as Eq. (7)
(೘)

(೙) మ

(೘) మ

(௫
ି௫ ) ା(௬
)
1
ି
ଶఙ మ
݁
ଶ
2ߨߪ
‫( ܩ‬௠,௡) =(௫ (೘) ି௫ (೙) )మ ା(௬ (೘) )మ
1
ି
ଶఙ మ
∑௖௫
௡ ୀ ଵ 2ߨߪ ଶ ݁(7)where cx is the number of nodes on the surface, x(m) and y(m) denote the x-coordinate and they-coordinate of node m. Then the shape identification introducing the Gaussian filter that parameter
σ is set 20 is carried out. The identified shape of cavity is shown in Fig. 8. As a result, the shape ofthe cavity becomes smoother and performance function at final iteration decreases in comparison
with the computational result without smoothing process (See Fig. 9). On the other hand, it is known
that the traction method is a smoothing technique that the displacement filed is obtained as velocity
field by using the gradient vector, elastic equation and boundary conditions[5]. Then elastic
coefficient E and Poisson ratio ν are given as 1 and 0, and gradient vector is smoothed from solutionof the elastic deformation analysis that the computational model is assumed as elastic body.
Computational results applying the traction method are shown below. Fig.10 shows the identified
cavity shape. It can be seen that there is no big difference in comparison with the identified shapewithout smoothing. In addition, performance function at final iteration is more than the case withoutsmoothing method (See Fig. 11).
4

Transactions on GIGAKU 4(1)(2017)04001/1-8

Fig. 6 Identified cavity without smoothing inCase A

Fig. 7 Variation of normalized performance
function without smoothing in Case A

Fig. 8 Identified cavity with Gaussian filter inCase A

Fig. 9 Variation of normalized performance
function with Gaussian filter in Case A

Fig. 10 Identified cavity with traction method
in Case A

Fig. 11 Variation of normalized performance
function with traction method in Case A

Moreover, shape identification analysis applying computational condition “Case B” are also
carried out. Figs.12 and 13 show the computational results without smoothing method. As a result,undulation of the inner surface can be seen in the identified cavity. Figs.14 and 15 show the cavity
shape at final iteration in computation applying the Gaussian filter and the variation of thenormalized performance function, respectively. It can be seen that the surface of the cavity has adimple similar to target shape. On the other hand, identification accuracy is not higher than the case5

Transactions on GIGAKU 4(1)(2017)04001/1-8

without smoothing method by performance function. And computational results using the traction
method are shown in Figs. 16 and 17. It is seen that the identified cavity surface is moreconcavo-convex than the result applying Gaussian filter (See Fig. 16). In addition, the value of theperformance function at final iteration decreases in comparison with the result without smoothing.

Fig. 12 Identified cavity without smoothing inCase B

Fig. 13 Variation of normalized performance
function without smoothing in Case B

Fig. 14 Identified cavity with Gaussian filter inCase B

Fig. 15 Variation of normalized performance
function with Gaussian filter in Case B

Fig. 16 Identified cavity with traction method
in Case B

Fig. 17 Variation of normalized performance
function with traction method in Case B

Finally, variations of normalized performance function for each conditions are shown in Figs. 18
and 19. In addition, Figs 20 and 21 show the comparison of the temperature at the observation point
6

Transactions on GIGAKU 4(1)(2017)04001/1-8

in each condition. In these figures, the temperature in each case is result at the final iteration. In CaseA, performance function at final iteration in iterative computation with the Gaussian filter is smallest
value. On the other hand, it is found that the performance function at final iteration can be
minimized by applying traction method in Case B. Similarly, in Figs 20 and 21, it is seen thattemperature at final iteration with the Gaussian filter is the closest to the observed temperature in allcases in Case A, and the temperature at final iteration with the traction method is the closest to theobserved temperature in all cases in Case B. Therefore, it can be said that the suitability of theapplication of each smoothing method depends on the shape of the cavity. Moreover, it appears that
if the appropriate smoothing method is employed in the cavity shape identification, the value of
performance function minimize, i.e., the temperature at the target point is close to the observed
temperature, in comparison with that without smoothing method.

Fig. 18 Comparison of normalized performance
function in Case A

Fig. 19 Comparison of normalized performance
function in Case B

Fig. 20 Time history of the temperature on theidentified shape in Case A

Fig. 21 Time history of the temperature on theidentified shape in Case B

5. Conclusions
In this study, the comparison of the smoothing technique employed in three-dimensional shapeidentification analysis for cavity in the resin structure was carried out. As for computational
condition, target shape was given as two type of model, and the Gaussian filter and the traction
method were employed as smoothing method. Results in this study are shown below.
(1) The Gaussian filter is effective for minimization of performance function in case of the target
shape with smoother surface.
(2) If the target shape is given as including a cavity with convex surface, the result applying thetraction method is better than that applying the Gaussian filter.

7

Transactions on GIGAKU 4(1)(2017)04001/1-8

Acknowledgments
This work was supported by the Program for High Reliable Material Design and
Manufacturingat Nagaoka University of Technology and a 2015 Grant for Joint Research between
Nagaoka College and Nagaoka University of Technology. The computations were mainly carried out
using the Fujitsu PRIMERGY CX400 computer facilities at Kyushu University’s Research Institutefor Information Technology.
References[1] T. Kurahashi and H. Oshita, Problem of Partial Reinforcement Corrosion in Concrete Using
Observed Temperature on Surface, Journal of The Japan Society for Industrial andApplied Mathematics, 21, 3 (2011) 221-240.
[2] H. Azegami and Z. C. Domain, Optimization Analysis in Linear Elastic Problems (Approach
Using Traction Method), Journal of The Japan Society of Mechanical Engineers, 60, 578(1994) 144-150.
[3] T. Kurahashi, K. Maruoka and T. Iyama, Numerical shape identification of cavity in three
dimensions based on thermal non-destructive testing data, Engineering Optimization, 1-17
(2016).
[4] H. Okuda and K. Deguchi, Gray Level Edge Behavior for the Gaussian Blurring Filters and Its
Application to Edge Extractions, Journal of Information Processing Society of Japan, 10,
36 (1995) 2244-2251.
[5] D. Fujii and J. Kakita, Shape optimization of frame structures using traction method, Journal of
Structural and Construction Engineering, 682, 77 (2012) 1881-1886.

(Received: 5 September 2016, Accepted: 4 April 2017)8

Transactions on GIGAKU 4(1)(2017)04002/1-7

Optimizing boundary surface of ceramics die in deep
drawing process
Hanh C. Nguyen1, Shigeru Nagasawa1,*, Kensei Kaneko2
1)

Nagaoka University of Technology, Nagaoka, Japan 940-2188
National Institute of Technology, Nagaoka College, Nagaoka, Japan 940-8532*E-mail: n_hanh@stn.nagaokaut.ac.jp

2)

The objective of the present study was to propose a new approach to optimize the geometry
parameters of ceramics die in deep drawing process. The optimization problem has been formulatedwith the objective of optimizing the maximum principal stress in the full stroke of deep drawing
process. In this study, FEM software was used to simulate the stress distribution on the die surface
and predict the brittle cracking occurrence on the die surface. Accordingly, the optimum ellipsoid
curvature on the curved surface of die hole can be obtained, and thus process design optimized.
1. Introduction

Due to its versatility, good strength and light weight design, deep drawing process is widely used invarious applications, especially in automotive and aerospace industries. However, after a time working, avariation of surface roughness, wear, and surface cracks gradually appear on the die surface relating to manyfactors: die radius, punch radius, friction between the die and the worksheet, and so on.
In a forming process of U-cup drawing, for an example as shown in Fig.1, a ceramics (Si3N4) die wasrecently tried to be used for pressing a 2mm thickness mild carbon steel disk, the diameter of which has
36mm, by using a 24.83mm diameter high carbon steel punch. In this case, due to a surplus high pressure on
the entrance surface of die, a fracture occurred at the die body after several times repetitions. Figure 2 showed
an example of cracked zone on the die surface. Observing the die surface with digital scanning microscope,
we can detect many cracking positions at the arc-end of ceramics die with Lsrf-dm varying from 8mm to 9mm as
was shown in the Fig.5. According to the stress analysis of worksheet without cracking on the die body, theMises stress of worksheet and the principal stress of die block were numerically simulated and the stress statesof them were relatively estimated. However, they can not predict the real fracture behavior neither its safety
designing condition.
upper punchworksheet

Ceramics die

HolderLower punch
Fig. 1 Deep drawing process model for making U-cup Fig. 2 A crack area on the ceramics die surface

Therefore, in this work, we propose a possibility and a methodology of prediction of cracking positon andhow to evaluate the modification of tool design. In order to develop a prototype FEM model composed ofelastic ceramic die body and a mild carbon steel in the axial symmetric condition, and also to predict the peak
position trajectory of the maximum principal stress and maximum shear stress, a case study of the real
production of U-cup shaped SPCE disk plate has been analyzed. This simulation model proposes an analysis
1Transactions on GIGAKU 4(1)(2017)04002/1-7

method for detecting the most severe point on a die surface, under specifying the radius of die hole and thecritical strength of the dies body. At the result, some optimization parameters can be found for designing theboundary surface of the ceramics dies.

2. Finite element simulation setup
The deformation behavior of cylindrical deep drawing process is simulated by using MSC. Marc version
2015 system - the Finite Element Method (FEM) software. This simulation model consisted of an elastoplastic body as a 2mm thickness worksheet of SPCE (mild carbon steel) and an elastic body as a ceramics die
block. The mechanical properties of SPCE worksheet was shown in Table 1. The in-plane mechanicalproperties of the SPCE worksheet were experimentally evaluated by the uni-axial tensile testing based on the
JIS-K7127. Figure 3 shows the relationship between the stress and strain under the tensile testing. The diebody of ceramics (Si3N4) was assumed to have a little plastic deformation at the breaking strength of 360MPa
which was gotten from the referenced tensile strength [1]. Generally, although the ceramics body such as
Si3N4 appears to be an elastic-brittle material, when a sort of damage criteria is diagnosed in the tool body, acertain plastic yielding is virtually necessary for evaluating the damage level of the tool body. Therefore, asmall breaking strain of eB= 0.0015 was virtually assumed here. This is used for detecting a severe-stressed
point from the aspect of damage function such as Oyane’s formula, Ayada’s formula [2].
Table 1 Mechanical properties

Worksheet
Die

Material

Young’s modulus
GPa

Poison’s ratio

Yield stress
MPa

Carbon steel
Si3N4

205
290

0.3
0.28

190

True stress σ, MPa

200

σY=190

True stress σ, MPa

σB=289300

σ=371.75e0.1783

100
eY=0.026000.2% offset

eB=0.283

0.1
0.2
True strain, e %Breaking stress, σB

400
350
300
250
200
150
100
50
0

σB=σY

0.05%
Breaking strain, eB
0

0.3

0.1

0.2

0.3

0.4

0.5

True strain e, %Fig. 3 Stress-strain curve of SPCE worksheet

Fig. 4 Breaking behavior of Si3N4

Figure 5 shows schematics of initial worksheet stacked on the upper level of die and the geometrical
relationship between the tool and the worksheet. Initially, the outside of working disk is attached to the entry
round surface of die and the upper punch is indented to the top side of working disk until the punch moves tothe stroke of 15.5mm. After releasing the upper punch from this pressing, the lower punch moves upward withthe same velocity as that of upper punch elevation.
Figure 6 shows the velocity diagram of upper and lower punches. An FEM simulation of the U-cup
drawing process was achieved in two stages. In the first stage, the upper punch moves down and acts upon the
carbon steel disk with a punching force. In the second stage, the upper punch was released, and the lower
punch moved up and acted upon the material with a punching force. In these two steps, both the upper punch
and the lower punch moved at a quasi-static state. The time-dependent thermal effect was not considered here.
For the sake of convenient execution, a virtual velocity of 2.5mm/s was assumed in MARC.

2

Transactions on GIGAKU 4(1)(2017)04002/1-7

18

2

Punch stroke, mm

12.415

Lsrf-dm

14.165

15
10
5
0
-5
-10
-15
-20
-25
0

10
15
Time, sUpper punch
Lower punch

Fig. 5 Schematics of worksheet geometries (unit: mm)and cracking positions in the experimentFig. 6 Punch motion diagram

13.65

5

3. Stress distribution on the ceramics surface
Figure 7 shows the variation of the peak maximum value of maximum1st principal stress σp1 (MPa) and
its occurred position Lsrf-p1 (mm) with respect to the upper punch displacement, while Fig.8 shows thevariation of the maximum shear stress t (MPa) and its occurred position Lsrf-ps (mm) with respect to theupper punch displacement. In this result, it is found that the peak of maximum shear stress occurred at the arcend inside position of round surface: Lsrf-ps = 8.82mm (on the un-loading process), and also the peak of firstprincipal stress occurred at the arc-end of round surface: Lsrf-p1= 10mm (on the loading process).
As a result, we can detect the critical point on the die body. It is useful for detecting the position of crack
initiation on the ceramics die surface. Comparing this result with the experimental cracked position (Fig.5),
we found that the cracked position was almost same as that of the peak maximum shear stress of t. The
position of σp1 was a little different position. Based on the FEM simulation, this work presented an analyzing
method for determining the peaks of the maximum principal stress distribution as well as the maximum shear
stress distribution on the die surface.
d/ts=13, Lsrf-p1= 10.0, σp1 = 337

16

350

14

300
d/ts=13250

12
Lsrf-p1

10

200

8

σp1

150

6

100

4
d/ts =15.550
ts=2mm

0
5

7
9
11
13
Normalized depth of indentation d/ts

015

Fig. 7 Peak of maximum 1st principal stress

3

2

17

Lsrf-p1 ,mm

Maximum principal stress σp1 ,MPa400

Transactions on GIGAKU 4(1)(2017)04002/1-7

Fig. 8 Peak of the maximum shear stressPrediction of severe position using a Damage model
Seeing the experimental result (Fig.2), cracks appeared on the die surface which could be because of manyfactors: die radius, punch radius, lubrication and so on. We used a method of predicting damage of drawing
die which is assumed to be applicable for ductile fracture. In this work, the Cockroft-Latham Damage
Indicator was considered in MARC [2] . It is often used in linear elastic fracture mechanics to predict brittle
fracture [3, 4]. The objective of this work is to provide a method for predicting the damage of drawing die
enabling to design of improved die surface by predicting brittle fracture under various conditions, andenabling us to design improved die surfaces. Fracture occurs when the damage factor has reached its criticalvalue of C. Here, we assumed that C=0.32 in the model.
஢eത̇ dt > C
∫ ౣ౗౮
ഥ
஢

(1)

Fig. 9 Damage occurs at the arc-end position of round surface

In the case of assuming eB=0.15% for die body breaking strain, the crack occurred at the arc-end inside
position of round surface: Lsrf-dm=9.4mm, d/ts=14.4., (on the un-loading process). The crack initiation
characteristics obtained via simulation are close to those got through the physical experiment as was shown in
Fig.5. The cracking position predicted by the simulation was shown in the Fig.9. This damage indicator isbased on the effect of Maximum principal stress over the plastic strain path. Therefore, when we define
eB=0.1% of breaking strain (perfect elastic linear behavior), this damage indicator can not detect the crack
point on the die surface. However, this damage model is applicable to predicting the damage for anotherfragile material such as GFRP laminated sheet and the acrylic sheet.
4Transactions on GIGAKU 4(1)(2017)04002/1-7

4. Parameter analysis of entry radius of die surface
Many factors affect the sheet metal forming process such as die radius, friction coefficient, blankholder force, and so on. Optimizing these parameters is an important task to reduce product cost. The dieradius has the greatest influence on the deformation behaviour of the deep drawing process (with an 89.2%
influence), followed by friction force with a 6.3% influence and the blank holder force with a 4.5% [5]. To
reduce the stress concentration on the die surface, we use appropriate ellipsoid curvature on the curved surfaceof die hole.

18

18

12,415
R6

12,415
R6

Rd

2

2

In this work, the authors presents a new approach to optimize the geometry parameters of die hole.
The optimization problem has been formulated with the objective of optimizing the maximum principal stressin the full stroke of deep drawing process. The current die surface (Fig.10 (a)) is provided an ellipsoid
curvature with rd set at 6mm, and Rd varying between 5mm and 8mm. The geometry is created in SolidWorks– a CAD/CAM software [6] and imported to the MARC system as showed in the Fig.10 (b).

R6

14,165

14,165a)b)rda. the current die surface b. variable ellipsoid curvature
Fig.10 Using appropriately curved surface of die holeFig. 11 Maximum principal stress and maximum shear stress in various die radius

5Transactions on GIGAKU 4(1)(2017)04002/1-7

To investigate the effect of entry radius on the stress distribution of die surface, different simulations
with different die radius were performed. Whereas the material properties, the punch profile, and the clearance
parameters are kept at the same conditions. The critical point occurs at the arc end of the die surface when thedie radius is varied from 5 mm to 8 mm. From Fig.11, it is found that maximum principal stress and maximum
shear stress of the drawn cups is decreasing as the die radius is increasing. However, when the die radius isincreased, wrinkle formation will occur in the worksheet [7]. Therefore, 6 mm in Rd is recommended as an
optimum value for this die fillet radius.
Punch force
Figure 12 shows the effect of the variation in die radius on punch force during the forming process,when the friction coefficent between worksheet and ceramics die is considered as Coulomb friction law withµ=0.23. When we increase the die radius, punch force tends to decrease. After the worksheet passes throughthe arc end of the die, the punch force suddenly increases because of the larger contact between the die surface
and the worksheet.

Fig. 12 Punch force required in various die radius

5. Conclusion(i) In order to diagnose the strength state of die body, the authors proposed to watch a state point of the peakmaximum principal stress and peak position (σp1, Lsrf-p1) in the die body for the full stroke motion. Plotting
this state point (σp1, Lsrf-p1) with respect to the punch stroke, the most severe position and stress value were
easily detected. It is found that this monitoring method is also useful for detecting an optimum design
parameters.
(ii) Comparing the simulated stress severe position and the experimental failure position, they were wellmatched. This results verified a validity of developed simulation model for discussing the strength of die bodymade by ceramics.
(iii) This study also proposes a method of predicting brittle damage for ceramics die in deep drawing process
and thus contributes to the design and parametric optimization of the ceramic dies.
Acknowledgments
This work was supported by a fund for developing a core of excellence as innovation and branding
project, from the GIGAKU Innovation Promotion Center, NUT, 2014-2015. The authors thank for providing acase of prototype ceramics die from Shimizu Industrial Co.Ltd.

6

Transactions on GIGAKU 4(1)(2017)04002/1-7

Nomenclature
Symbolσp1 or σmaxt

Lsrf (p1,ps,dm)
d/ts
࣌
ഥ
തė
C,
R
Rd, rdDefinition
Maximum principal stress, MPa
Maximum shear stress, MPa
Arc length (with 1st principal, shear maximum, damage), mm
Normalized depth of indentation, mm
Effective von Mises stress, MPa
Effective elastic-plastic strain rate
Material threshold for damage
Radius of circle of entry hole in die
Radius of elipsoid, in the long axis and short axis

References[1]
Ferro-Ceramic
Grinding,
I.,
Wakefield,
Mass.,"Material
Properties
Charts,"https://www.ceramicindustry.com/ext/resources/pdfs/2013-CCD-Material-Charts.pdf.
[2] MSC Software, Marc® 2010 Volume A: Theory and User Information.
[3] Bouchard P.O., B. F., Chenot J.L., "Modeling of sheet metal cutting by coupling damage and crack
propagation theories," Proc. simulation of Material processing: Theory, Methods and Applications, pp. 10011006.
[4] Goijaerts A. M., G. L. E., Baaijens F. P. T., 1999, "Prediction of Ductile Fracture in Metal Blanking,"
Journal of Manufacturing Science and Engineering, 122(3), pp. 476-483.
[5] Padmanabhan, R., Oliveira, M. C., Alves, J. L., and Menezes, L. F., 2007, "Influence of process parameterson the deep drawing of stainless steel," Finite Elements in Analysis and Design, 43(14), pp. 1062-1067.
[6] MSC Software, Release Overview - Marc 2014.
[7] Kakandikar G.M., N. V. M., 2007, "Optimization of Forming Load and Variables in Deep Drawing Process
for Automotive Cup Using Genetic Algorithm," On-line Journal, http://www.optimizationonline.org/DB_HTML/2007/03/1606.html.
(Received: 7 October 2016, Accepted: 6 April 2017)7

Transactions on GIGAKU 4(1)(2017)04003/1-9

Bending behaviour of polypropylene sheet subjected to two-line wedge indentation

Akihiro Yamamoto1*, Satawat Singprayoon2, Shigeru Nagasawa1, Seksan Chaijit3
Dept. of Mechanical Engineering, Nagaoka University of Technology, Nagaoka, Japan
2)
Exchange research student in NUT from Pathumwan Institute of Technology, Bangkok, Thailand
3)
Dept. of Mechatronics Engineering, Pathumwan Institute of Technology, Bangkok, Thailand*E-mail: akiyamamech@gmail.com
1)

In this paper, the bending characteristics of a 0.3mm thickness Polypropylene (PP) sheet
subjected to a sharp two-line wedge indentation were investigated. Through the experiments, the formability of the scored zone was characterized with the groove depth andfolding (normal and reverse wise) direction, and compared with that of Polyethylene
Terephthalate (PET).
1. Introduction
Half-cutting of worksheets is a well known process and used for producing packaging boxes, especially for folding thick resin sheets [1,2,3]. The relationship between the mechanical conditions forthe half-cutting and the folding profile of the worksheet are important characteristics to create a successful crease in a packaging box. The process of wedge indentation of a worksheet has been formally explained by Hesse [4] and Mizuguchi et al. [5] in their respective works on technology for
packaging sheets. The most important mechanical conditions regarding the half-cutting are the cutting wedge profile, the indentation depth of the wedge and the material properties of the worksheet.
On the other hand, the folding characteristics of wedged (scored) worksheet are mainly influenced
by the deformation profile of the scored zone. This folding characteristics can be verified by studying the bending moment response regarding a certain folding angle. Empirically, the residual thickness after scoring (thickness of the worksheet minus the indentation depth of the wedge cutter) of theresin sheet is important in determining the quality of folded profile and keeping a certain strength on
the hinge left after folding [3]. However, these folding characteristics seem be sensitive and changewith respect to the material properties. Thus, there is a need to experimentally reveal the bending
behavior of a worksheet with respect to the wedge indentation depth and the folding direction (normal and reverse wise).
Nagasawa et al. [6,7] have reported about the effects of scored residual thickness on the bending
characteristics of a 0.5mm thickness Polyethylene Terephthalate (PET) board subjected to a two-line
wedge indentation. Here, a laminated structure of three layers (G-PET/A-PET/G-PET) was investigated and the strength of laminated interface was discussed with respect to the scored residual
thickness. The gradient of bending moment resistance for folding and the in-plane tensile strength of
the grooved zone were also discussed for the specified PET sheet. However, this approach was not
used to investigate other resin sheets.
Recently, thick sheets of Polypropylene (PP) are becoming more popular for producing transparent carton boxes. Although there are reports of usage and fundamental properties [1,8], and re-1

Transactions on GIGAKU 4(1)(2017)04003/1-9

search articles have been reported for cutting strength [9,10], there are not almost any academic reports regarding the half-cut folding of PP sheets. Therefore in this work, Polypropylene SG-170 (PP)sheets of 0.3mm thickness were chosen to investigate the formability and the bending characteristicsof the scored groove in respect to the folding direction and residual indentation depth.
2. Experimental method
Figure 1 shows schematics of experimental apparatus for scoring on a 0.3mm thickness PP
sheet. In the experimental apparatus, a facet blade (two-line wedge) mounted on the upper crosshead
moved downward with the feed velocity V=0.03 mms-1. The two-line wedge was composed of thefirst (apex) angle α=90° and the secondary α’=30°. Here, a height of the first facet h was 0.15 mm.
The wedge blade was made of cemented carbide (FM10K) and it had a tip thickness of 0.8 μm,
height of 9 mm, length of 30 mm, thickness of 0.9 mm and the surface roughness Ra=0.09 μm (as agrinding surface).

(b) Details of blade profile and definition of d,dR, h, α, α’Fig. 1 Schematics of experimental apparatus for scoring(a) Layout of experimental apparatus

In general, since there is an equivalent-spring effect on the blade displacement, the indentation
depth d can be principally estimated with the difference between the blade displacement and device
elongation. In this work, the device elongation was not inspected. Alternatively, the indentation
depth of wedge blade d= t-dR was calculated from the measured residual thickness dR after half cutting (scoring). The applied pushing force F was measured by a load cell.
The PP sheet was a rectangular specimen with a length of 50 mm and width of 15 mm. The position of the blade was vertical to the worksheet, while the angle f of the cutting line direction withrespect to Machine direction (MD) of the worksheet was chosen as f =90°.
Table 1 In-plane mechanical properties of PP
Young’s Modulus
Yield StressTensile Strength
Breaking Strain
E / GPaσY / MPaσB / MPaεB / 2.09
23.26
222.401.76
Before indenting the blade to, the worksheet was gently placed on the counter plate without any
fixture. The depth of wedge indentation was considered as a half cutting condition: d/t= 1-dR/t =0.5~0.95. Regarding a shallow indentation of the wedge blade, when 0<d/t<0.4, due to a low bending
stiffness of PP sheet and its large tensile strength, the folding of scored zone behaved unsuccessfully.
Therefore, the shallow indentation of wedge blade was not considered in this work. Here, the folding test
was carried out using a bending strength tester CST-J1 [11] explained in the below. Table 1 shows the2Transactions on GIGAKU 4(1)(2017)04003/1-9

in-plane mechanical properties of the PP sheet. Here, the strain rate of feeding was chosen as 0.1 s-1.
Comparing the mechanical properties between the PP sheet and the reported PET sheet [6], the ratioof the Young’s modulus, yield stress, tensile strength and breaking strain were as follows: the ratio of
E.PP/EPET = 0.92, the ratio of sY.PP/sY.PET = 0.33, the ratio of sB.PP/sB.PET =1.47 and the ratio of

eB.PP/eB.PET = 1.48. Namely, the PP sheet is easily flown out but difficult to be broken out, comparedto the PET sheet.
The scored PP sheet was bent in either the normal or reverse direction using a crease strength
testing machine CST-J1 [11] for measuring the bending moment per unit width M N and the rotation
angle θ° as shown in Fig. 2. The distance between the load cell and the scored position (the center ofrotation) was set for L=10mm. The scoring and folding experiments were carried out under the following conditions: room temperature of 297 K, room humidity of 50 %RH.

(a) Normal bending(b) Reverse bending

Fig. 2 Schematics of folding test of scored sheet
The work sheet was bent in either the normal or reverse direction with rotation velocity of
ω=0.1 rps (revolution per second). Here, the rotation angle θ was repeatedly (twice) varied from zero
up to the right angle (0~90°). Figure 3 (a) and (b) illustrate representative response models of relationship between the bending moment resistance M and the rotation angle θ for the normal wise andreverse wise, respectively.
0.5

d/t=1-dR/t=0.7

0.4
0.3
0.2

q1

0.1

Mp1
0
-10

qp1 q2

Bending moment resistance
M /Nm m-1

Bending moment resistance
M /Nm m-1

0.5
H2

qc

H1

C1

C2

10
30
50
70
Rotation angle q /°90

d/t=1-dR/t=0.7
0.4
0.3

H2

0.2
0.1

qp1

q1

Mp1
0
-10

H1

q2
C1

C2

10
30
50
70
Rotation angle q /°90(a) Normal bending(b) Reverse bending
Fig. 3 Response model of bending moment diagram and parameters for characterizing its response.
Through the bending moment response, the first break or inflection bending moment Mp1, itscorresponded angle θp1, the first gradient C1, the second gradient C2, the initial starting position (angle) θ1, the second starting position (angle) θ2, and the two kinds of work hardening gradient H1, H2
were measured when varying the indentation depth 1-dR/t for the normal and the reverse wise bend3

Transactions on GIGAKU 4(1)(2017)04003/1-9

ing. The number of samples was five pieces for each case. Here, in the case of normal wise bending,the gradient H1 was estimated as the average gradient between the first break point Mp1 and thestarting point of interference, while the gradient H2 was estimated as that between the middle pointof 75° and the tracking (the maximum) position of q=90°. In the case of reverse wise bending, empirically choosing the middle position as q=60°, the gradient H1 was estimated between the Mp1 andthe middle point of 60°, while the gradient H2 was estimated between the middle point of 75° and thetracking position of q=90°. In the case of normal wise folding, since there was an interference ofwedged groove, the interference based bent inflection angle qc was recorded.

Fig. 4 Sectional views of half-cut (scored) zone of PP sheet(a) 1-dR/t=0.50, normal(b) 1-dR/t=0.80, normal(c) 1-dR/t=0.50, reverse(d) 1-dR/t=0.80, reverse
Fig. 5 Side views of scored PP sheet when folding up to θ = 90 °3. Experimental results and discussion
Figure 4 shows photographs of sectional views after scoring with 1-dR/t= 0.50~0.80. The residual (permanent) wedged angle A was measured from those photographs. The value of A was discussed later with the interference resistance during the folding test. Figure 5 shows some photo4

Transactions on GIGAKU 4(1)(2017)04003/1-9

graphs of the scored PP sheets when the worksheets were bent in either the normal or reverse direction with the right angle q = 90°. Depending on the scored depth, the profile of scored part (groove)
was remarkably altered. Figure 6 showed the first break moment Mp1 with respect to the normal andreverse bending direction. The relationship between Mp1 and 1-dR/t and that of qp1 were linearly approximated as Eq.(1)-Eq.(4) for the normal and reverse directions.
Mp1 = 0.41-0.47 (1-dR/t) (for normal)(1), Mp1 = 0.40-0.47 (1-dR/t) (for reverse) (2)qp1= 26.5-22.0 (1-dR/t)(3),(for normal)

qp1= 27.2-24.9 (1-dR/t) (for reverse) (4)Comparing the first break (or peak maximum) of Mp1 of PET [11], Mp1 of the PP sheet andAngle of first break bending
moment qp1 deg.

First break of bending moment
Mp1 Nm m-1

Eq.(1),(2) were less than 0.1 times of that of the PET sheet, while qp1 of the PP sheet and Eq.(3),(4)were about a half of that of the PET sheet.

0.25
0.20

Eq.(1) normal

0.15
0.10

Eq.(2)
0.05 reverse
0.0
0.4
0.6
0.8
1.0
Blade indentation depth 1-dR/t(a) First peak bending moment

25.0
20.0

Eq.(3) normal

15.0
10.0

Eq.(4)
5.0 reverse
0.0
0.4
0.6
0.8
1.0
Blade indentation depth 1-dR/t(b) Angle at first peak bending moment

Fig. 6 Comparison of first peak point with respect to normal and reverse bending
According to the theory of plastically collapse moment with a beam bending [12], the collapse(ultimate) moment per unit width Mpc is estimated as Mpc = sY∙t2/4 = 0.53 Nm∙m-1 when 1-dR/t=0.

The nominal collapse (ultimate) moment Mpc = sY∙dR2/4 was estimated as 66~70 % of experimental

Mp1 when 1-dR/t=0.5~0.8. This result is a little small but seems to have a primary contribution. Using the yielding stress ratio sY.PP/sY.PET = 0.33 and the thickness ratio 0.3/0.5= 0.6, the value of Mpc
of 0.3mm PP sheet is estimated as 0.11 times of that of 0.5mm PET sheet. It was found that the orderof collapse moment matched to the range of the first break (or peak maximum) bending moment
when varying the indentation depth of 1-dR/t. Here, since the ratio of Young’s modulus of PP by
PET was E.PP/EPET = 0.92, the strain for starting the yield seems to have the same order as that of

sY.PP/sY.PET=0.33. Since the ratio of q1p for both cases was about 0.5, the yield ratio of 0.33 was notequal to this break angle ratio, but they were related with each other.
Figure 7 shows a representative relationship between M and θ when varying the indentation
depth of blade. The interference of scored groove occurred in the normal wise bending. As the result,a remarkable increase of bending moment M was detected. Through the results of relationship between M and θ, the gradients of C1, H1 , H2 and C2 were arranged with 1-dR/t. The first starting angle q1 was almost zero for all the cases in this work.
Figure 8 shows the first gradient C1 and second gradient C2 of bending moment resistance withrespect to the normal and reverse bending direction. In the case of the normal folding, C1 and C2
5

Transactions on GIGAKU 4(1)(2017)04003/1-9

were linearly approximated with 1-dR/t as Eq.(5), Eq.(6). Here, the relation of 0<C2<C1 seems to becaused by damages of folding repetition. It was found that C1 of the normal bending was almost
same as that of the reverse bending, and also C2 was similar between the normal and reverse bending.
The ratio of C1, C2 between PP and PET [6] was about 0.13~0.19 at c/t=0.5. Since the first gradient

0.7
0.6
0.5

0.7

1-dR/t=0.50
0.60
0.70
0.80
0.90

Bending moment resistance M (Nm)m-1

Bending moment resistance M (Nm)m-1

C1 is approximately expressed as Mp1/qp1, the values of C1, C2 seem to be strongly affected by theyield stress of the worksheet and the residual thickness dR.

0.95

0.4
0.3
0.2
0.1
0
-20

0

20

40

60

1-dR/t=0.50
0.60
0.70
0.8
0.90
0.95

0.6
0.5
0.4
0.3
0.2
0.1
0

80 100-20

Rotation angle, θ °

0

20

40

60

80 100

Rotation angle, θ °Gradient of bending moment
C1, C2 Nm m-1deg-1(a) Normal bending(b) Reverse bending
Fig. 7 Relationship between bending moment resistance and rotation angle when varying indentation
depth of blade.

0.016
C1, normal
Eq.(5)
0.014
C1, reverse
0.012
C2, normal
0.010
0.008
C2, reverse
0.006
0.004 Eq.(6)
0.002
0.000
0.4
0.6
0.8
1.0
Blade indentation depth 1-dR/t

Fig. 8 Relationship between first and second gradients of bending moment and indentation depth ofblade, along with normal and reverse bending modes.
C1= 0.029-0.03(1-dR/t)≈0.03dR/t (for normal)(5)

C2= 0.021-0.023(1-dR/t)≈0.02dR/t (for normal)(6)

Fig. 9 (a) shows the first half gradient of bending moment H1, while Fig. 9 (b) shows the second
half gradient of bending moment H2 with the normal and reverse bending. The difference of the latter gradient H2 between the normal and reverse bending were obviously caused by the groove interference of normal bending. There is a small difference in the H1 for 1-dR/t≤0.6. It seems to becaused by the interference behavior in the groove. In order to discuss this problem, the bent inflec6

Transactions on GIGAKU 4(1)(2017)04003/1-9

tion angle qc and wedged angle A were shown in Fig.10. Here, the plotted data were evaluated fromthe average and the representative photograph. The value of qc was almost equal to that of A, namely
the initial clearance of groove was sufficiently kept when 1-dR/t≥0.7. In case of shallow blade indentation, the wedged groove clearance is insufficient and the bent inflection of qc becomes invisible.
When the groove depth is the deeper, the difference of H1 between the normal and reverse bending
becomes elastic deformation without interference and similar with each other. In the case of PET
sheet, the wedged angle A was relatively larger than that of PP sheet. For an example, A= 84 degrees

0.012
0.010
0.008
0.006
0.004
0.002
0.000
0.4

Gradient of bending moment
H2 Nm m-1deg-1

Gradient of bending moment
H1 Nm m-1deg-1

when 1-dR/t=0.57 [6].

normalreverse

0.6

0.8

1.0

Blade indentation depth 1-dR/t

0.012
0.010
0.008
0.006
0.004
0.002
0.000
0.4

normalreverse

0.6

0.8

1.0

Blade indentation depth 1-dR/t(a) First half gradient H1(b) Second half gradient H2

60

Second starting position (angle)
q2 / deg.

Angle of wedged and interferenceA, qc / deg.

Fig. 9 Anaphase characteristics of bending moment
80

Residual,wedged angle A

PET, reverse

PET: Ref.[6]

60

40
20

PET, normal

40

Bent inflection angle qc

20

0
0.4

0.6
0.8
Blade indentation depth 1-dR/t

1.0

Fig.10 Wedged angle and inflected-bent position

PP, reverse
00.4

PP, normal

0.6
0.8
1.0
Blade indentation depth 1-dR/t

Fig.11 Second starting position (angle)

Regarding the repeated (as twice folding) elastic behavior, the second gradient C2 and the second
starting position q2 were investigated. Figure 11 shows the comparison of q2 with the normal andreverse folding, along with that of PET [6]. Here, the graph was derived from the average for eachcase. Seeing Eq.(5),(6), C2 was about 67% of C1. This relation was almost the same level with thenormal and reverse folding. Seeing Fig.11, it was found that q2 was remarkably affected by the folding direction (normal and/or reverse) in the case of PP sheet for 1-dR/t≤0.8. When choosing
1-dR/t≥0.9 in the case of PP sheet, the value of q2 was 60-70 degrees (67~78% of 90 degrees),
which was almost the same as that of PET sheet.
4. Conclusions
In order to detect allowance of mechanical conditions of two-line wedge blade, the crease characteristics of polypropylene sheets of t= 0.3 mm thickness were investigated. Regarding the normal
7

Transactions on GIGAKU 4(1)(2017)04003/1-9

and reverse bending of scored part, the bending characteristics such as the first/second stiffness gradient C1, C2, the first/second work hardening gradient H1, H2, the first break (or peak maximum)
bending moment Mp1, and its corresponded angle of θp1 were experimentally investigated with respect to the indentation depth of blade. Some of them were compared with that of PET sheet.
Through this experiment, the following results were revealed.
1) Bending stiffness parameters such as C1, C2, Mp1, qp1 were characterized with the yielding
stress sY and the residual thickness dR. They were strongly related to the plastic collapse moment.
2) The profile of wedged groove was different from that of PET sheet, when using the same tip
shape of two-line wedge blade as for 0.5< 1-dR/t <0.95. Depending on the wedged angle A, in thecase of the normal folding, a bent inflection point of bending moment occurred at the rotation angle

qc ≈ A for 1-dR/t ≥ 0.7.
3) Regarding the repetition (twice folding), the relationship between the second starting angle q2
and the residual thickness dR was different from that of PET sheet. It tends to be similar to that of
PET sheet for the range of 1-dR/t >0.9.
When choosing 0.5<1-dR/t <0.8, the tip profile of the wedge blade is ought to be changed and its
scoring characteristic should be investigated furthermore.
References[1] R.J. Crawford, Plastics Engineering, 3 eds., Burlington: Elsevier Butterworth-Heinemann,(1999), pp.1-40.
[2] Japan Society for Technology of Plasticity eds., in: Processing Technology of Plastics, Bending,(2016), pp.229-230 (in Japanese).
[3] Tadashi Hashimoto, Processing Technology of Transparent Packaging Box and Related Problems, CARTON-BOX, vol.23 (2004), pp.44-45 (in Japanese).
[4] F. Hesse, H.J. Tenzer, Grundlagen der Papier-verarbeitung, VEB Verlag fur Buch und Biblio-thekswesen, Leipzig, (1963), pp. 58-60 (in Germany).
[5] S. Mizuguchi, eds, Advanced Technology on Industrial Packaging and Transportation, Fuji
Techno System, Tokyo, (2002), pp. 58-59 (in Japanese).
[6] Shigeru Nagasawa, Kiminori Sato, Mitushiro Murayama, Yasushi Fukuzawa, Effect of Tool
Condition on Bending Characteristics of Polyethylene Terephthalate Board Subjected to Two-Line
Wedge Indentation, Proceedings of the 10th International conference on technology of plasticity(ICTP 2011), Aachen, Germany, pp.984-989.
[7] Shigeru Nagasawa, Cutting and Bending of Plastics Sheet, Journal of Japan Society for Technology of Plastics, vol.57 no.665 (2016), pp.867-872 (in Japanese).
[8] Hideaki Ishihara, Fundamental and Application of Plastics Films, Converting Technical Institute,Tokyo, (2010), pp.33-52 (in Japanese).
[9] Satoshi Okuda, Woo Sik Choi, Mechanism of Cutting Fracture of Polypropylene by Impact on
Wedge-shaped Target, Journal of Chemical Engineering of Japan, vol.14 no.2 (1981), pp.149-153.
[10] W. Bauer, D. Wustenberg, Fracture Behavior of Polypropylene under Dynamic Cutting andShearing Action in Granulators, Chemical Engineering and Technology, vol.25 no.11, (2002),
pp.1047-1051.
8

Transactions on GIGAKU 4(1)(2017)04003/1-9[11] Katayama Steel Rule Die Inc., “Crease Stress Tester CST-J-1”, (accessed on September, 2016),
http://diemex.com/sale/cst.html.
[12] Hideaki Kudo, Plasticity, Morikita publishing, Tokyo, (1968), pp.16-19 (in Japanese).
(Received: 10 October 2016, Accepted: 17 February 2017)9

Transactions on GIGAKU 4(1)(2017)04004/1-6

Improvement to Reduce Delay Time in Transportation
to Disaster Area using Similarity to Electric Circuit

1)

Chata Yorvarak1, a, Thanat Rungsirathana2, Kazumasa Takahashi1,
Toru Sasaki1, Takashi Kikuchi3, and Nob. Harada1

Dept. Electrical, Electronic and Information Technology, Nagaoka University of Technology,
2)
School of Engineering, University of the Thai Chamber of Commerce,3)
Dept. Nuclear System Safety Engineering, Nagaoka University of Technology,1603-1 Kamitomioka-machi, Nagaoka 940-2188, Japan*E-mail: yorvarak026@gmail.com, nob@nagaokaut.ac.jp

Transportation is closely related to many sectors such as economic, politicaland social environment of the country. It is necessary to identify suitabletransportation method for keep the pace of the country’s economy and life styleof citizens in disaster times. Currently most frequent disaster in the world isflooding and this study focused to reduce the delay time in transportation for
those crisis to sending disaster reliefs. Similarity between the logistics andelectrical circuits will be used to stimulate the transportation models which canbe used in disaster times. Simulation work in this study was carried in thisstudy using the MATLAB/Simulink software and parameters were introduced
to system considering the similarity between logistics and electrical circuits.
Initially study focused on water distribution to disaster areas in chaos time
where parameters such as amount of water required were identify as amount ofcharge in electrical circuits for the simulation. Also requirement of number oftucks for the delivery was identified as number of pulses in electrical circuitand capacity of water storages were identified as capacitance of the electricalsystem, etc. Centralized distribution model was stimulating for understanding
the behaviors of the electrical charge flow and change of parameters which aresimilar to electrical circuit shows changes in the flow of charges with time inthe electrical circuits which can be adopted to logistic patterns. Present study
results show decentralized patterns is more suitable than the centralized
patterns for improve to reduce the delay time in transportation in disaster areas.

1. Introduction
A disaster incident is defined as incident of any natural or man-made, that result in levels ofserious damages to human, properties, economic or environmental losses which affecting the people
directly [1]. Disaster occurrence have increased significantly in last several years in many countries, inthat flood disaster has the highest occurrence frequency among all natural disasters. It is necessary tohave a plan to flow model to help decide how to rapidly supply aid to victims of a disaster within thiscontext [2]. Include find the optimization route to suitable in these situations because reservation
others alternative can find the optimize route to supply as soon as possible for respond people who are
in disaster. The main purpose is to provide a faster response, possibly by decreasing delay time
between the origin and destination of supplies. We propose an optimization model that could be usedwithin a decision-aid tool are currently responsible for providing rapid responses in disaster relief.
However, logistic can offer lots of logistic patterns and it is essential to identify and stated proper
logistic plan to send relief items to disaster areas. In these situations, disaster relief logistics (DRL)can be applied [2], [3].
When using DRL and developing logistic systems, it is necessary to understand the processes andpathways of that plan. Using electrical circuit can be helpful to run and tested the models before goesfor the actual logistic plan. Through the concepts of the electrical circuits, with the similarity withlogistics, can tested the different scenarios of the disaster and can identify the most suitable plan forthe cause. For testing logistic plans through electric circuits can be stimulate through different
software programs and it will reduce the collection of huge amount of data and testing time of the1

Transactions on GIGAKU 4(1)(2017)04004/1-6

logistic plans. MATLAB/Simulink is one of the powerful applications to stimulate electrical circuits
for optimization, statistics, mathematical functions, or data analysis, then there is a very quick learning
curve for using it in image processing. This study focused on optimization model for find the best way
to quickly respond to disaster victims because of the situation is often unclear in the immediate
aftermath of the disaster. The option of a route that can respond most and to identify the similarity
between the logistic and electrical circuits and stimulate developed electrical circuit models to identify
the suitable logistic patterns to assess the logistic models for disaster reliefs.

2. System Parameters and Modelling
Proper logistic network can reduce the delivery time and also the cost of distribution reliefs indisaster time. However, in this study electrical circuits will be analyzed using the MATLAB/Simulink
software and it requires to introduce proper parameters for analysis. Fig. 1 shows the simple concept
for logistic model to deliver water to disaster area point “1” from the distribution center/warehouse
which is point “0”. Introduced electrical circuit network for the logistic model is shown in the Fig. 2
and Table 1 shows the relation between the logistic and electrical circuit with set parameters forevaluation. For the evaluation model water can be considered to be electric charge.

Fig. 1 Basic logistic model from the point “0” to the destination point “1”.
Table 1 Similarity relation between logistic and electric circuit for the case of water supply.

Fig. 2 Logistic equivalent electrical circuit for the simple logistic model for water delivery model.
2Transactions on GIGAKU 4(1)(2017)04004/1-6

3. Assumptions-water supply/current source
Typical source of water delivery is shown in Fig. 4 where 6 trucks with each capacity of 4000L.
So it looks like rectangular pulse series as shown in Fig. 4.

Fig. 4 Assumed for source storage.
For the simulation water supply will be consider as a digital pulse with a delay time which was
calculated using the assign parameters of distance between two nodes and speed of truck. For the
analysis of the two scenarios which are 1. Storage without any initial items and 2. Storage with initialamount of items were considered

4. Analysis of Electrical Circuit using Simulink Model
4.1 Basic Model
The model for a logistics using electrical circuit was created by using MATLAB/Simulink software asshown in Fig. 5. The process of connection between the blocks-cells and flow between each block
created through Simulink. For analysis, two models were created where the Model 1 concidering no
initial storage and the Model 2 with 50% of initial storage. Simulink model was created using thecomponents of its library and parameters were defiend in the property boxes of each component.

Fig. 5 Basic logistic model created using MATLAB/Simulink.
3Transactions on GIGAKU 4(1)(2017)04004/1-6

4.2 Assigning Values for Simulation
Primary values were assigned as shown in the Table 2 for run the basic analysis of the two models
and relevant platform were used to assign values to the program for simulation. In this simulation
program one liter were considered to be one micro coulomb (1 l = 1 μC) and one hour to be one
second (1 h = 1 s). Typical simulation results for these models were shown in Fig. 6. We could see that
initial delay of water supply for the both cases but in black curve (Model 1) there is a crisis without
the water till receive its first batch of water. On the other hand, red curve (Model 2) shows that water
can be supplied from initial storage till receive its initial water supply. Rate of supplying water wasconsidered as the function of water in the tank and results shows with higher storage of water can help
to increase the supply rate for disaster victims.
Table 2 Assigned parameters
Parameter
Number of trucks

Model 18

Model 2
6

Truck capacity (million liters)
Distance (km)Truck speed (km/h)
Storage of tank (million liters)

0.025
450
65
0.6

0.025
600
65
0.6

Initial storage (million liters)

0.0

0.05

Fig. 6 Simulation results for the Model 1 and Model 2.
4.3 Changes in parameters4.3.1 Same set of parameters values for both models
Fig.7. shows that the frequency of both curves were same and the delay time also same as the
input parameters are similar to each other. However, for the result 2 which had initial storage (redcurve) always went over the result 1 (black curve) as its storage facility contained more relief items
every time, comparing to both cases. The distribution rates somewhat different compare to each otheras red curve to black curve as the result 2 had more initial storage.
4.3.2 Change in distance between point 0 and point 1
Change in the parameter of distance between distribution point and disaster point showed how it
affected the delay time for distribute relief from point 0 to point 1. Longer distance in the result 2
compare to that in the result 1 caused larger delay time which was shown in Fig. 8. As otherparameters were the same, frequency should be the same between two nodes.
4Transactions on GIGAKU 4(1)(2017)04004/1-6

Fig. 7 Evaluation for same parameters (black curve - result 1, red curve - result 2).

Fig. 8 Evaluation for change in distance (black curve - result 1, red curve - result 2).

Fig. 9 Evaluation for change in number of delivery trucks.
(black curve - result 1, red curve - result 2)5Transactions on GIGAKU 4(1)(2017)04004/1-6

4.3.3 Change in number of delivery trucks
Influence of change in number of trucks on the delivery curves was shown in Fig.9 as no storage
situation overcome the initial storage as the number of trucks bringing more capacity of the disaster
reliefs per a day in the result 1 compare to that in the result 2. As other factors were kept constantfrequency dependence could be important point in disaster area which need higher priority.

Fig. 10 Proposed model for reduce delay time in disaster relief transportation.
Conclusion
Logistic related electric circuits are relatively easy for analyzing through the simulation than thereal time data collection and analysis. This study was carried out with MATLAB/Simulink simulate
software as it provides easy platform for analysis and it helps to make variation for the possible realtime scenarios. For the preliminary results, it was observed that the amount of water supply per a day
had an impact as how much relief could be sent at disaster time. When sending the disaster relief items,
longer the distance increased the delay time of delivery, and therefore, initial capacity tank and/or
initial storage of the disaster areas played a major role in addressing the disasters as soon as possible.
So, in this study we were recommended to have another sub-terminal between main distribution points
to disaster area as show in Fig.10. and built some storage or tank to keep initial storage for distribute
relief items till the enough relief items coming from distribution centers.in order to reduce delay time.
both the closer distribution point and the larger initial storage could increase the more delivery ofdisaster relief items.
References[1] Mei-Shiang Chang, Ya-Ling Tseng, Jing-Wen Chen, "A Scenarioplanning approach fo theflod emergency logistics preparation problemm under uncertainty," Transportation Research ,vol. 43, no. E, pp. 737-754, 2007.
[2] A. Thomas & M. Mizushima, "Logistic Training: necessity or luxury?. Forced Migration
Review," Fitz Institute, 2005.
[3] Mohammad Rezaei-Malek, Reza Tavakkoli-Moghaddam, Behzad Zahin, Ali Bozorgi-Amiri,"An interactive approach for designing a robust disaster relef logistics network with perishable
commodities," Computers & Industrial Engineering, vol. 94, pp. 201-215, 2016.
(Received: 7 October 2016, Accepted: 2 December 2016)6

Transactions on GIGAKU 4(1)(2017)04005/1-7

An Observation on Potential Alkali Silica Reactivity of Natural
Aggregate in Myanmar
Pyae Phyo Kyaw1, Thynn Thynn Htut2,*
Quality Control Department, Tetlann Concrete Co., Ltd.,
2)
Building Department, Taisei Myanmar Co., Ltd.
No. 32, The Tokyo Enterprise Building, 6 ½ Mile, Hlaing Tsp, Yangon, Myanmar*E-mail: thynnhtut@taiseimyanmar.com
1)

Most aggregate used in concrete are almost chemically inactive, however,mineral constituents of some aggregates react with alkali hydroxide, which
originates mainly from the Portland cement. The chemical reaction, internallyoccurred alkali–silica reaction (ASR) in either concrete or mortar under acertain condition, e.g. high humidity and sea saline water, may result
deleterious expansion of concrete or mortar. The prior evaluation ofalkali-aggregate reactivity by the reliable test is required for the utilization ofnatural aggregate in concrete for assurance on the longer service life of thestructures.
The objective of this study is to identify the potential alkali–silica reactivityof several aggregates which able to buy in Yangon region, Myanmar. To
achieve this objective, four numbers of fine aggregate quarry which areoperated to produce river sand was selected and investigated. In addition,
three numbers of coarse aggregate quarries were also chosen to investigate
potential alkali–silica reactivity.
Four samples of river sand, one sample of river shingle and three samples ofcrush stone were collected from the mentioned quarries and the accelerated
Bar Mortar test was conducted on these samples to assess the potential
alkali–silica reactivity. According to ASTM C 1260 suggested method, three
Bar Mortar tests were performed on each sample. Based on the obtained
results, it was found that four samples of all seven samples were considered
as non-reactive due to their average expansion rates which are less than 0.1%.
Two samples which were collected from Wan Be Inn and Ohn Chaw quarry
are considered as moderately reactive which average expansion rate between
0.1% and 0.2%. Only one sample taken from Kamaryut quarry is concluded
as reactive with the average expansion rate of 0.3429% which is more than
0.2%.

1. Introduction
Generally, alkali–silica reaction (ASR) is considered as an internal chemical reaction between
the alkaline components in the cement and active silica–silicate based mineral constituents of someaggregates. The use of aggregates susceptible of ASR in concrete can cause severe degradation of
the concrete structures because of several reasons e.g. high humidity and sea saline water [1].
Map-like cracks appear when concrete is suffered from ASR [2]. Many researchers have stated that
ASR is generated due to combine action of existence of a significant quantity of reactive silica inaggregates, having critical level of alkalis (sodium and potassium) in pore solution of concrete andwater or moisture from the external source. Utilization of potentially alkali–silica reactive aggregatesneeds to conduct reliable tests for evaluating their alkali–silica reactivity [3,4].
Generally, it is considered that aggregates from the quarries in Yangon area are not reactive.
The portion of aggregate quarry which may have potential reactive aggregates need to be identified.
Hence, the research study undertaken has following aims and objectives:· To identify the potential alkali–silica reactivity of several aggregates available in Yangon
region, Myanmar
1Transactions on GIGAKU 4(1)(2017)04005/1-7· To identify Coarse aggregate which meets criteria of AASHTO.
To achieve this objective, four numbers of the fine aggregate quarry which is operated to
produce river sand was selected to investigate alkali–silica reactivity of sand. In addition, threenumbers of the coarse aggregate quarry were also chosen for the same purpose. Aggregate samples
are collected from the mentioned quarry and an extensive laboratory tests: accelerated Bar Mortar
test was conducted on these samples to assess the potential alkali–silica reactivity.

2. Evaluation of Alkali-Silica Reactivity
2.1 Fine and Coarse Aggregate Quarry
The selected quarries for evaluation of alkali-silica reactivity of aggregate are listed inTable 1. The aggregate manufactured from those quarries are mainly supplied to the constructionindustry in Yangon region.
Table 1 List of sand and coarse aggregate quarry
Sr.
No

Quarry Name

Material

Name of River

RegionTownship1

Wartayar

River sandHlaing River

YangonShwe Pyi Thar

2

Thilawa

River sandYangon River(Tidal Mouth)

YangonThanlyin3

Kamayut

River sandYangon River

YangonKamayut

4

Wan Be Inn

River sand

Sittaung River

Mon state

Thein Za Yat

5

Mi Kyaung YeRiver shingle

Ayeyarwaddy

Bago

Pyay

6

Ohn Chaw

Crush stone-

Mandalay

Ohn Chaw

7

Moke Pa Lin

Crush stone-

Mon state

Moke Pa Lin

2.2 Sample Collection

Thilawa

Kamayut

Wan Be Inn

Ohn Chaw

2

Moke Pa LinTransactions on GIGAKU 4(1)(2017)04005/1-7

Fig 1 Aggregate Samples
The aggregate samples are collected from the quarry stated in Table 1 and those samples areshown in Fig 1. The coarse aggregate is crushed and the crushed aggregate is poured onto the sieves
to separate them into their various sizes. Then it is washed each size with a water spray over thesieve to remove adhering dust and fine particles from the Aggregate. Each portion is placed into theoven at 105 Degree Celsius and allows to dry overnight. Each such portion is individually stored in aclean container provided with a tight-fitting cover. Just before casting, the aggregate is combined in
the following grading requirements as in table 2 for aggregate with relative density ≥2.45 g/cm3.
Table 2 Grain Size Distribution for Mortar Bar [5]Sieve Size

Mass (%)

Mass (g)2.36 mm10991.18 mm

25

247.5

600 µm25

247.5

300 µm25

247.5

150 µm15

148.5

Passing
4.75mmRetained on

2.36mm
1.18mm
600µm
300µm

2.3 Moulding Mortar Bar
Portland cement that meets the requirements of ASTM Specification C 150 is used. Lumps are
removed and cement are passed through an 850 -µm sieve before use. Then, aggregate and cement
are mixed in a mixing bowl. The moulds are filled up with two approximately equal layers, eachlayer being compacted with the tamper. The mortar is placed into the corners, around with gauge
studs and along the surface of the moulds with the tamper until a homogeneous specimen is
obtained.
2.4 NaOH Solution Preparation
Initially, using a weighing balance with weigh 40g of NaOH pellets and dissolved in 900ml ofdistilled water. Then, it pours the solution into a 1000ml volumetric flask and fills up the solution to
the point with more distilled water.
2.5 Conditioning
The temperature of moulding room and dry materials are maintained between 20 °C and27.5 °C. Also, the moist room and closet are controlled to be 23°C±1.7°C. The relative humidity ofthe moulding room is kept at 50%. The specimens are stored in the storage oven that controlledtemperature at 80 ± 2 °C.

3Transactions on GIGAKU 4(1)(2017)04005/1-7

Fig 2 Mortar Bar Mould

Fig 3 Storage Oven

2.6 Initial Storage and Reading
Each mould is placed in the moist cabinet immediately after moulds have been filled. The
specimens are kept in the moulds for 24 ± 2 h. The specimens are removed from the moulds and,while they are being protected from loss of moisture, give the name to specimens and make an initialcomparative reading. The initial and all subsequent readings to the nearest 0.002 mm are made andrecorded. The specimens are placed in a storage container with sufficient tap water to totally
immerse them. The container is sealed and placed in an oven at 80.0 ± 2.0°C (176 ± 3.6°F) for aperiod of 24 hours. The container is removed from the oven then, the specimens are taken off fromthe container and the surfaces of mortar bar are dried with a towel paying particular attention to thetwo metal gage studs. The zero reading of each bar is made immediately after drying and read assoon as the bar is in position. As specified in ASTM C-1260, the process of drying and reading iscompleted within 15 ± 5s of removing the specimen from the water. After comparative readings on
the remainder of the bars, all specimens made with each aggregate sample are placed in a container
with sufficient NaOH, at 80.0 ± 2.0°C (176 ± 3.6°F) for the samples to be totally immersed. The
container is sealed and returned to the oven. The subsequent comparator readings of the specimens
are made periodically, with at least three intermediate readings, for 14 days after the zero reading, atapproximately the same time each day.
2.7 Calculation of Percent Expansion
The difference between the zero-comparative reading of the specimen and the reading at eachperiod to the nearest 0.001 % of the effective gage length is calculated and record as the expansionof the specimen for that period. The average expansion of the three specimens of a given
cement-aggregate combination to the nearest 0.01 % as the expansion for the combination for agiven period is recorded. The expansion (En) of each mortar bar from the following equations.
௅ ି௅‫ܧ‬௡ = ೙ଶହ଴ మ × 100

and, if required, the value of thermal expansion‫ܮ‬ଶ − ‫ܮ‬ଵܶ‫= ܧ‬
× 100
250
En = expansion at the time n days, after immersion in the NaOH solution at 80°C in percent [%]TE = thermal expansion after heating to 80°C in the water bath [%]Ln = comparator reading at time n days after immersion of mortar bars in the NaOH solution [mm]L1 = Initial comparator reading [mm]L2 = comparator reading after heating to 80°C in the water bath [mm]3. Results and Discussion
3.1 Aggregate Reactivity Classification
Aggregate reactivity is classified based on the 16 days expansion values in accordance with4

Transactions on GIGAKU 4(1)(2017)04005/1-7

Table 3, which tabulates the requirements of ASTM standard specification.
Table 3 Aggregate reactivity classification [5]Mortar Bar Expansion (%) in NaOH (80%)

Classification

Less than 1%

Non-reactive

Less than 2%

Reactive

More than 2%Moderate reactive

Individual results shall not differ from the mean by more than 15%.
3.2 Experimental Results
Fig. 3 shows expansion of mortar bar incorporated with fine aggregate from War Ta Yar
Quarry. The expansion curves for WTA-1, WTA-2 and WTA-3 are in the same trend. The expansion
rates are very small. All specimens have under 0.03 % in the expansion of 16 days.
0.4

0.4

0.35

0.35
0.3
WTA-1

0.2

WTA-2

0.15

WTA-3

Expansion [%]Expansion [%]

0.3
0.25

TLW-2

0.2

TLW-3

0.15

0.1

0.1

0.05

0.05

0

TLW-1

0.25

0
0

2

4

6

8

10

12

14

16

18

0

2

4

6

Age of Mortar Bar [days]

8

10

12

14

16

18

Age of Mortar Bar [days]

Fig 3 Expansion of mortar bar (Wartayar)

Fig 4 Expansion of mortar bar (Thilawa)

In figure 4, the expansion of mortar bars composed of sand from Thilawar quarry shows small
expansion in early of the experiment. Afterward, gently increase in expansion happened started from
6 days till the age of 16 days. The expansion of all specimen at the age 16 days don’t reach 1%.
Expansion curves of KMT-1, KMT-2 and KMT-3: sand from Kamaryut quarry
incorporated into mortar bars are shown in Fig.5. The sudden increase in expansion found between
10 and 12 days of the experiment. It can be seen in Fig 10 that the average expansion curve for
Kamaryut sand pass through a moderately reactive zone and it reached reactive zone at age of 12days. Then, expansion rates significantly increase until the end of the experiment: the last test for theage of 16 days. It is clear that sand from Kamaryut quarry is the reactive aggregate.
Expansion of specimen by Wan Be Inn aggregate is shown in Fig. 6. Expansion curves for all of
the mortar bars: WBL-1, WBL-2 and WBL-3 gradually increase from the beginning of theexperiment to test on 12 days. After that moment, the expansion rate doesn’t increase and become
stable until the end of the experiment.

0.4

0.4

0.35

0.35
0.3

Expansion [%]Expansion [%]

0.3
0.25
0.2
0.15

KMT-1

0.1
0.05

0.25
0.2
0.15

WBI-1

KMT-2

0.1

WBI-2

KMT-3

0.05

WBI-3

0

0
0

2

4

6

8

10

12

14

16

0

18

2

4

6

8

10

12

14

16

18

Age of Mortar Bar [days]

Age of Mortar Bar [days]

Fig 5 Expansion of mortar bar (Kamayut)

Fig 6 Expansion of mortar bar (Wan Be Inn)
5

Transactions on GIGAKU 4(1)(2017)04005/1-7

Fig 7 Expansion of mortar bar (Mi Kyaung Ye)

Fig 8 Expansion of mortar bar (Ohn Chaw)

Expansion of mortar bar made with river shingle from Mi Kyaung Ye quarry is shown in Fig 7.
Expansion % at the age of 6 days slightly grow up for all mortar bar reach 0.03% expansion. At theage 10 and 12, all specimen exhibit expansion and MKY-1 and MKY-3 reach 0.035% expansion.
Afterward, MKY-02 has the great increase in expansion at the age of 14 days and it subsequently
goes stable at the end of the experiment: 16 days in age of mortar bar. The tendency of expansion
curves for all: MKY-1, MKY-2 and MKY-3 are the same. It can be said that the results are reliable.
Expansion of mortar bar made up of aggregate from Ohn Chaw Quarry is displayed in Fig.8. In
its experiment, all mortar bar: OHC-1, OHC-2, and OHC-3 are in the same tendency in expansion
from the beginning of experiment till the end of the experiment.
As shown in Fig. 9, expansion % of mortar made up of sand from Moke Pa Lin quarry increase
with a great extent until the age of 8 days counted from the start of the experiment. The highest
expansion % of MPL-3 is 0.045% at the age of 16 days and the others are 0.04% at the same age.

0.4
0.35

Expansion [%]

0.3
MPL-1

0.25

MPL-2

0.2

MPL-3

0.15
0.1
0.05
0
0

2

4

6

8

10

12

14

16

18

Age of Mortar Bar [days]

Fig 10 Compared Mortar Bar Expansion and
ASR Reactivity Range

Fig 9 Expansion of mortar bar (Moke Pa Lin)3.3 Compared Mortar Bar Expansion and ASR Reactivity Range
It can be clearly seen in Fig.10 that the average expansion rate for Wartayar sand throughout theexperiment is very low. So that, fine aggregate from War Ta Yar quarry can be deduced that
non-reactive aggregate.
An average expansion of Thilawa sand doesn’t reach ASTM’s limit of the moderately reactive
zone until the end of the experiment. So that, sand from Thilawar sand can be said as non-reactive
aggregate.
In Fig 10, the average expansion curve for Kamaryut sand pass through the moderately reactive
zone and it reached reactive zone at age of 12 days. Then, expansion rate significantly increases until
the end of the experiment: the last test for the age of 16 days. It is clear that sand from Kamaryut
6

Transactions on GIGAKU 4(1)(2017)04005/1-7

quarry is reactive aggregate.
The average expansion of Wan Be Inn sand at 16 days is 0.184% and that is very close to reach
ASTM’s reactive zone. It can be said that aggregate from Wan Be Inn quarry is moderately reactive
based on the current result.
Average expansion % of MKY-1, 2 and 3 are presented in Fig.10 incorporate with aggregate
reactivity classification in accordance with ASTM C-1260. ASTM C-1260 classified that 0 to 1% of
expansion of mortar bar of 16 days is non-reactive. So that it can be deduced that river shingle from
Mi Kyaung Ye quarry is non-reactive.
The specimen made of aggregate from Ohn Chaw quarry showed the significant increase in
expansion from the start to 10 days and those reach 0.1% expansion at the age of 10 days. The
average curve of expansion rate stays lower region of the moderately reactive zone in Fig. 5.4. So
that, it is concluded that crush aggregate from Ohn Chaw quarry is moderately reactive.
For the mortar bar made up of Moke Pa Lin aggregate, it can be concluded that aggregate form
Moke Pa Lin quarry is non-reactive due to the average of expansion for MPL-1, 2 and 3 shown in
Fig. 10 stay at the middle of ASTM C-1260 non-reactive zone.

4. Conclusion
In order to identify potential alkali-silica reactivity of fine and coarse aggregate, seven numbersof aggregate quarries located in Myanmar was selected and studied. In the study area, three types ofaggregate were observed in which main rock type is naturally found in the study area and thefollowing fact can be concluded.
(a) River shingle from Mi Kyaung Ye quarry is non-reactive aggregate.
(b) Crush aggregate from Ohn Chaw quarry is moderately reactive aggregate.
(c) Fine aggregate from War Ta Yar quarry can be deduced that non-reactive aggregate.
(d) Fine aggregate form Moke Pa Lin quarry is non-reactive aggregate.
(e) Sand from Kamaryut quarry is reactive aggregate.
(f) Sand from Thilawar sand can be said as non-reactive aggregate.
(g) Aggregate from Wan Be Inn quarry is moderately reactive.
References[1]. M. Berra et al., “Alkali-silica reactivity criteria for concrete aggregates,” Materials andStructures, vol. 38, no. 3, pp. 373–380, 2005.
[2]. J. A. Farny and B. Kerkhoff, Diagnosis and control of alkali-aggregate reactions in concrete.
Skokie, IL: Portland Cement Association, 2007.
[3]. F. Rajabipour et al., “Alkali–silica reaction: Current understanding of the reaction
mechanisms and the knowledge gaps,” Cement and Concrete Research, vol. 76, pp. 130–
146, Oct. 2015.
[4]. Z. Owsiak, “Alkali–aggregate reaction in concrete containing high-alkali cement and granite
aggregate,” Cement and Concrete Research, vol. 34, no. 1, pp. 7–11, Jan. 2004.
[5]. ASTM standard C 1260, 2014, “Standard Test Method for Potential Alkali Reactivity ofAggregates (Mortar-Bar Method)” ASTM International, West Conshohocken, PA, 2014,
DOI: 10.1520/C1260-14, www.astm.org.
[6]. U.S.Departemnt of Transportation, Federal Highway Amimistration, M. D. A. Thomas et al.,
Selecting Measures to Prevent Deleterious Alkali -Silica Reaction in Concrete: Rationale
for the AASHTO PP65 Prescriptive Approach.

(Received: 8 October 2016, Accepted: 5 January 2017)7

Transactions on GIGAKU 4(1)(2017)04006/1-8

Parallel implementation of Entropic lattice Boltzmann method for flow
past a circular cylinder at high Reynolds number
Ayurzana Badarch1,*, KhenmedekhLochin2, Hosoyamada Tokuzo3
1)
Graduate School of Engineering, Nagaoka University of Technology,
2)
School of Applied Science, Mongolian University of Science and Technology, Mongolia
3)
Department of Civil and Environmental Engineering, Nagaoka University of Technology
1603-1 Kamitomioka-machi, Nagaoka 940-2188, Japan*E-mail: ayur_426@yahoo.com

Lattice Boltzmann method (LBM) has been receiving enormous attention
from researcher to solve fluid related phenomena. LBM in fluid dynamics
simulation, one meet some difficulties, for instance it becomes unstablefor high Reynolds number flows and requires computational time andmemory for large scale simulations. To avoid this, we used here theEntropic LBM (ELBM) of Karlin’s group [7] and implemented parallelcode on graphical processing units (GPU). For parallel computation, weused GPGPU system of Nagaoka University of technology, which isequipped with Tesla M2050 processors. To verify accuracy and stability,
we have solved double periodic shear layer flow in serial computation.
The accuracy and stability of simulation by ELBM were superior to thesimulation of standard LBM. Simulations of flow past a circular cylinder
is carried out by parallel ELBM, where Reynolds number varies up to
140000. Using the GPU in simulation, computation time speeds up until
10 times faster than that of using central processing units (CPU). The
results show that the parallel code of ELBM can solve two-dimensional
turbulent flow in arbitrary geometry at the higher stability condition
without using any other stabilization techniques.

1. Introduction
The simulation of traditional computational fluid dynamics deals with Navier-Stokes equation (NSE).
Unlike this, lattice Boltzmann method (LBM) solves the discrete Boltzmann kinetic equation in its
simplified form, known as Bhatnagar-Gross-Krook (BGK) equation. Founded in the eighties of lastcentury, LBM become recently a widely used method of computational fluid dynamics (CFD). In
LBM, the fluid is simulated by molecular velocity distribution functions with discrete velocities onregular lattice. There are two operations on the distribution functions: streaming and collision. The
fluid molecules propagation is described by streaming, which is non-local operation against the local
collision operator, describing the collision between molecules [1]. LBM has its intrinsic feature: easy
to parallelize on computer [2]. On the other hand, one of long standing problem in fluid flow
simulation is the simulation of turbulent flow at high Reynolds number. In CFD, the turbulent flow
simulations are usually performed either by computationally expensive direct numerical simulationof NSE or by solving averaged NSE [3, 4, 5, 6]. Direct application of the LBM for the turbulent flow
simulation leads to instability of computation, because of the low viscosity of fluid. Laminar flows at
low Reynolds number (Re) become turbulent when Re reached the value of Re=300 in the case offlow past a bluff body. One of the promising method of stabilizing the LBM turbulent computation isthe entropic LBM (ELBM) of Karlin’s group, in which, the distribution functions are straightened tosatisfy the maximum condition of the entropy at the every time step of simulation. In this paper, first
we confirm accuracy and performance of the ELBM solving double periodic shear layer flow. Thenwe perform the calculations of high Reynolds number turbulent flow past a circular cylinder using
the ELBM of Karlin et al [7, 8] in parallel algorithm.

1

Transactions on GIGAKU 4(1)(2017)04006/1-8

2. Numerical model
Recently, a simple entropic stabilizer was proposed to stabilize LBM simulation at highReynolds number in incompressible flow. The entropy of the lattice Boltzmann equation (LBE) canbe written as଼݂௜ܵ[݂] = − ෍ ݂௜ ݈݊ ൬ ൰,
‫ݓ‬௜

(1)

௜ୀ଴

where ݂௜ is the density distribution functions, ‫ݓ‬௜ is the weight and ܵ[݂] is the entropy. Further wewrite the ݂௜ as a component of the velocity moments:
݂௜ = ݇௜ + ‫ݏ‬௜ + ℎ௜ ,(2)where ݇௜ is the kinematic part of populations, which depends only on the locally conserved fields,‫ݏ‬௜ is the shear part of populations, which depends on the stress tensor and ℎ௜ is the higher ordermoments of populations as a linear combination of the remaining higher order moments. Introducing௘௤௘௤
deviations ∆‫ݏ‬௜ = ‫ݏ‬௜ − ‫ݏ‬௜ and ∆ℎ௜ = ℎ௜ − ℎ௜ , the extremum condition of the ܵ[݂] take the form:଼෍ ∆ℎ௜ ln ቆ1 +௜ୀ଴(1 − ߚߛ)∆ℎ௜ − (2ߚ − 1)∆‫ݏ‬௜
݂௜

௘௤

ଵ

ቇ = 0,
ଵ(3)

where ߚ can be defined from the kinematic viscosity as ߥ = ܿ௦ଶ ቀଶఉ − ቁ and the maximum ofstabilizer ߛ can be computed asଶ଼

1
1 〈∆‫∆|ݏ‬ℎ〉
ܺ௜ ܻ௜
ߛ = − ൬2 − ൰with 〈ܺ|ܻ〉 = ෍ ௘௤ .
ߚߚ 〈∆ℎ|∆ℎ〉
݂௜
∗

௜ୀ଴

(4)

Non equilibrium part of the distribution function can be written as௘௤݂௜ − ݂௜ = ∆‫ݏ‬௜ + ∆ℎ௜ ,(5)and it is used for calculation of higher order moments ∆ℎ௜ . The deviations of shear part of the
populations are defined by following velocity moments:‫ݑ‬ଶ − ‫ ݒ‬ଶ෩ ௫௬ − ቈ෩ቇ ,⎧∆‫ݏ‬଴ = ߩ ቆ4‫ݒݑ‬П቉ܰ
2⎪⎪
ߩ 1 + ܿ‫ݔ‬௜ ‫ ݑ‬+ ‫ݑ‬ଶ − ‫ ݒ‬ଶ෩ − [2ܿ‫ݔ‬௜ ‫ ݒ‬+ 4‫]ݒݑ‬П෩ ௫௬ ቇ ,቉ܰ∆‫ݏ‬ଵ,ଷ = ቆቈ
⎪
2
2
(6)
ߩ −1 − ܿ‫ݕ‬௜ ‫ ݒ‬+ ‫ݑ‬ଶ − ‫ ݒ‬ଶ⎨෩ − [2ܿ‫ݕ‬௜ ‫ ݑ‬+ 4‫]ݒݑ‬П෩ ௫௬ ቇ , and∆‫ݏ‬ଶ,ସ = ቆቈ
቉ܰ
⎪
2
2
⎪
−‫ݑ‬ଶ + ‫ ݒ‬ଶ − ܿ‫ݔ‬௜ ‫ ݑ‬+ ܿ‫ݕ‬௜ ‫ݒ‬
ߩ
⎪∆‫ݏ‬
෩
෩ቇ ,
቉ܰ
ହ,଺,଻,଼ = ቆ[4‫ ݒݑ‬+ ܿ‫ݔ‬௜ ܿ‫ݕ‬௜ + 2ܿ‫ݕ‬௜ ‫]ݑ‬П௫௬ + ቈ
42⎩
෩ are the central moments and ܿ‫ݔ‬௜ and ܿ‫ݕ‬௜ are the discrete velocities of lattice.
෩ ௫௬ and ܰ
where ПThe central moments are defined by the natural moments and velocities as
଼

଼

଼

௜ୀ଴

௜ୀ଴

௜ୀ଴

෩ ௫,௬ = ෍ ݂௜ ܿ‫ݔ‬௜ ܿ‫ݕ‬௜ − ‫ ݒݑ‬and ܰ෩ = ෍ ݂௜ ܿ‫ݔ‬௜ଶ − ෍ ݂௜ ܿ‫ݕ‬௜ଶ − (‫ݑ‬ଶ − ‫ ݒ‬ଶ ).
П

(7)After defining stabilizer with above formulas, the post collision state of distribution function iscomputed in the form:
݂௜ᇱ = ݂௜ − ߚ(2∆‫ݏ‬௜ + ߛ ∗ ∆ℎ௜ )(8)and the streaming is done by copying distribution functions to their designated directions. The
equilibrium distribution functions are defined by the Maxwell distribution as‫܋‬௜ ∙ ‫܋( ܝ‬௜ ∙ ‫)ܝ‬ଶ ‫ܝ ∙ ܝ‬௘௤݂௜ = ‫ݓ‬௜ ߩ ቈ1 + ଶ +
− ଶ቉,
(9)
ܿ௦
2ܿ௦ସ
2ܿ௦
where ߩ is the density, ‫܋‬௜ is the discrete velocity and ‫ ܝ‬is the macroscopic velocity. The densityand velocity are estimated by the moments of the distribution functions asேߩ = ෍ ݂௜ ,
௜ୀ଴

2

ே

ߩ࢛ = ෍ ࢉ௜ ݂௜ .
௜ୀ଴

(10)

Transactions on GIGAKU 4(1)(2017)04006/1-8

The collision expression given in Eq. 8 is called KBC operator in Entropic LBM [7], while standard
LBM uses BGK operator for collision.

3. Parallel implementation of the Entropic LBM
In this study, standard and entropic LBM are implemented in serial and parallel way using
Fortran 90 programming language. A stable computation of the ELBM is gained by paying time andmemory as the computational cost, which will be explained in the Section 5. At the other points,
engineering problems are often characterized as a large scale or long term phenomena in space andtime. LBM solution for these problems in CPU based simulation ended up at difficult challenges
because of the computational cost. The implementation of the ELBM on GPU can solve the two
important difficulties of computation: instability and high computational cost.
Nagaoka University of Technology (NUT) has a GPGPU system, which is available for students
and staffs for their research activities. The GPGPU system is powered by 16 nVIDIA Tesla M2050
graphic processing units, which has ability to perform 562GFLOPs in double precision. Also CUDA
Fortran version 11.3 has been installed on the GPGPU system. Grasping this technical potential andnumerical method, we have made an attempt to implement efficient parallel code for fluid flow onGPU. In CUDA programming, code has to be devoted to two particular sections named host codeand device code [9], as shown in Fig. 1. Host code refers that code executes on the host CPU
machine, whereas device code refers that code executes in the GPU machine. The codes use
memories on dedicated machines, i.e., variables used in device code must use a GPU memory.

Fig. 1 General architectural scheme of the implemented parallel code on the GPGPU system.
As the main program in the CUDA FORTRAN code is executed on the host, all variables in both
host and device codes need to be declared specific memories of their use in the host code.
Macroscopic variables, distribution functions and population parts are declared on the global
memory (device memory) and could only be used in device code. As purpose of the ELBM code,
initial condition and loading geometric data are implemented on the host code, in which device andhost codes need to exchange data. So called kernels, subroutines of Fortran language, implemented
for device code is launched from the host code. On the right hand side of Fig. 1, the host code
composing the kernels with used equations is given, which is performed on the predetermined blocks
and threads in GPU architecture. An entropic part of the numerical algorithm tested as separated
kernels, for instance a kernel for population part and a kernel for stabilizer. The separated kernels use
low memory consumption, however, total computation time of two kernels were higher than theirunified version in one kernel. The defining population parts and stabilizer is optimized to use low
memory on the local memory of the GPU in one kernel, since the declaration of array in the kernel
use a lot of registration, which might be caused the abortion of the kernels. Besides the streaming
step, all required equations and formulas in the ELBM are easily implemented in parallel way. To
perform the streaming in efficient way, two parallel distribution functions are used and they3

Transactions on GIGAKU 4(1)(2017)04006/1-8

synchronize in the host code after the calculation of one time step. To improve performance ofparallel code, some optimization of array declaration as ݂(‫ݔ‬, ‫ݕ‬, ݅) for distribution function insteadof the regular array structure of distribution function ݂(݅, ‫ݔ‬, ‫ )ݕ‬is applied. Storing the constantparameters used in both code in their memory save small amount of the computational time. A single
lattice is assumed to be a thread and a computational domain has divided as blocks horizontally.
4. Numerical test on doubly periodic shear flow
Doubly periodic shear layer flow is often considered a benchmark [7,10] case for under-resolved
simulation of smooth flows with sharp features. Shear flow had been solved by serial code ofstandard LBM and ELBM in resolutions of N=128, 256, 512, and 1024. Initial conditions for flow
field are given by‫ݕ‬⎧‫ݑ‬଴ ‫݊ܽݐ‬ℎ ቆ݇ ቀ − 0.25ቁቇ , ‫ܰ ≤ ݕ‬/2,⎧⎪ܰ
⎪‫= ݑ‬⎪‫ݕ‬⎨(11)⎪‫ݑ‬଴ ‫݊ܽݐ‬ℎ ቆ݇ ቀ0.75 − ቁቇ , ‫ܰ > ݕ‬/2,
ܰ
⎨
⎩
⎪
‫ݔ‬
⎪
‫ݑߜ = ݒ‬଴ sin ቆ2ߨ ቀ + 0.25ቁቇ ,ܰ⎩where ݇(= 80) is the parameter controlling the width of the shear layer and ߜ(= 0.05) is theparameter creating small perturbation of velocity in y-direction, which initiates a Kelvin-Helmholtz
instability. The turning over time of the shear layer is defined as ‫ݐ‬௖ = ܰ/‫ݑ‬଴ , where ‫ݑ‬଴ (= 0.04) is
the initial velocity defining Reynolds number for flow as ܴ݁ = ‫ݑ‬଴ ܰ/ߥ. In numerical test, Reynolds
number was set at 30000. The simulation results of vorticity at ‫ݐ‬௖ by two methods are shown in Fig.
2 in order of grid resolutions. The simulation of standard LBM was failed before the convergence of
the under-resolution at low resolution of domain, for instance unstable appears at t = 2200 before thevortex roll-ups in the interface of shear layers, shown in top figure of Fig. 2 (a).

Fig. 2 Double peroidic shear layer flow simulations on the different grid resolutions. Top figures
were produced by using the standard LBM, while bottom figures were produced by using theentropic LBM. Vortex roll-ups appears at (a) tc=3200, (b) tc=6400, (c) tc=12800 and (d) tc=25600.
As an increase of grid resolution, relaxation time of the standard LBM increases and the computation
become rather stable. The computation successfully survived on the further cases of simulations in
the case of standard LBM. However, as shown in top of Fig. 2 (b), two small additional roll-ups
created at linear part of vortex field, which is the source of unstable solutions. As expected shapes ofvortex, entropic LBM produces flow field at the even small resolution of grids and additional
4Transactions on GIGAKU 4(1)(2017)04006/1-8

roll-ups emerging from unstable of the numerical solution does not appear in any case of simulation.
This stable simulation can be performed at very small kinematic viscosity due to the coarse grid andhigh Reynolds number. The stabilizer in ELBM is a self adaptive local parameter and impacts on the
resolving of the flow fields. The distribution of the stabilizer and its value range do not depends on
the Reynolds number. For computational cost of the serial codes, ELBM code was 2 to 4 times
slower than the standard LBM code.

5. Parallel performance of the ELBM code
The parallel code for ELBM is implemented by the way that uses lower memory and efficientcomputation, as explained in the Section 3. We have solved flow past a cylinder in laminar regime(Re = 200) to evaluate performance of the parallel code, since the standard LBM can’t survive athigh Reynolds number. The simulation domain formed as 300 grids in length and 100 grids in width.
The simple bounce back boundary condition was imposed for top wall, bottom wall and a circularobstacle. The Zou/He boundary condition was applied with a characteristic velocity of ‫ݑ‬଴ =0.1 forinlet and the zero gradient boundary condition is used for outlet boundary. PGI Fortran compiler on
the GPGPU system is used to perform parallel and serial computations to evaluate parallel
performances. Time criteria to stop simulations was set at t = 50000 in all simulations.

Fig. 3 Percent of elapsed times used by the
kernels in single computational time of the
ELBM in double precision.

Fig. 4 Speedup ratio and parallel efficiency of
the ELBM code with regards to the number of
threads.

As investigating deeply the performance of the parallel ELBM code, Entropic part (labeled “Ent”
in Fig. 3), computing stabilizer and population parts, consumes around 41% of the all required timeto perform kernels, shown in Fig. 3. The second largest consumer was defining macroscopic
variables and synchronization of two groups of distribution functions (labeled “Macr” in Fig. 3),
since the streaming process done by using the two groups of distribution functions.
The computational domain was divided as blocks in horizontally and threads in a block. A
number of blocks can be chosen. Figure 4 shows speedup ratio and parallel efficiency of parallel
computation on GPGPU with respect to the number of threads used in simulation. The parallel
efficiency is determined as the ratio between the speedup ratios to the warp sizes of GPU. The
parallel code can speed up until 10 times faster than its serial version on the same machine. The
highest speedup ratio was found at the 64 threads in a block. The computational efficiency would be
convenient where the number of threads is lower than the maximum grid number of domain. In theother words, the number of threads should be smaller than the maximum number of grid inhorizontal or vertical direction.
Further we examined computational throughput and memory bandwidth for four cases of grid
resolutions, shown in Fig. 5. The number of thread was 512 in the simulations. Computational
throughput and memory bandwidth have linear relation with the grid resolutions. Comparing to the
maximum performance of the GPU, our code performs at the middle rate. Actually, memory band
width was inefficient, uses small amount of the possible memory bandwidth for the data transfer and5Transactions on GIGAKU 4(1)(2017)04006/1-8

management. Therefore, effective usage of the memory band width for parallel ELBM code needs tobe improved. Elapsed time against the increase of grid resolutions plotted in Fig. 6 shows
performances of the code on the different grid resolutions. Hence, the increase of the computationaltime has parabolic relation with the increases of grid resolution. Steepness of the parabola can be
decrease if we increase effective memory bandwidth in the code. To improve effectiveness of thecode, possible kernels can be merged in such a way that uses low memory in local memory of thedevice. Interestingly, LBM index, which is index to show parallel performance of the LBM in theway that how much lattice is updated per second (LUPs), increased as the grid resolution increases
in Fig. 6.

Fig. 5 Computational throughput and memory
bandwidth against the grid resolutions

Fig. 6 The elapsed time and efficiency of ELBM
computation against the grid resolutions

6. Numerical simulation of flow past a cylinder
Flow past a bluff body has been widely studied both numerically and experimentally [13] at highReynolds number [3-6]. The enhanced version of LBM using sub-grid scale [11, 12] and direct
numerical simulation like entropic approaches [7, 8] are the main contributions for those studies. We
have used a domain used in [12] for our simulation. The boundary conditions are the same explained
in section 5. Preliminary test simulations were performed for the Re=1000 and Re=3900 in order totest accuracy. We defined Strouhal number and compared with the other studies conducted in same
condition using different techniques in Table 1.
Table 1 Strouhal numbers of flow around circular cylinder at different Reynolds number
Sources/methods of the results

Strouhal number

Re = 1000 – Present study

0.217

Extended LBM (using Smagorinsky sub-grid scale model) [12]0.208

Experiment [13]0.212Standard k-epsilon model [5]0.148Re = 3900 – Present study

0.211

Extended LBM (using Smagorinsky sub-grid scale model) [12]0.215

Experiment [13]0.215Standard k-epsilon model [5]0.171Re = 140000 – Present study

0.228

Extended LBM (using Smagorinsky sub-grid scale model) [12]0.209

The Strouhal numbers defined by ELBM was in good agreement with the results of experiment
6

Transactions on GIGAKU 4(1)(2017)04006/1-8

and other simulation methods. In Fig.7, the vorticity field of fully developed 2D turbulent flow at
simulation time t = 10000 at Re = 140000 performed by parallel ELBM is visualized.

Fig. 7 Vorticity field of flow past cylinder at Re=140000 on coarse grid.
A performance of the ELBM can clearly describe Reynolds number effect on the cylinder. The
flow field can be captured exactly comparing with performance of LBM extended by Smagorinsky
model. It was revealed that Smagorinsky model for LBM overestimates eddy viscosity [14] at highReynolds number and it leads underestimates of velocity field. Unlike it, ELBM does not modify
viscosity in the simulation and corrects velocity field by the maximum entropic condition. For coarse
grid, ELBM can resolve flow field and generate vortex in sub-grid. Figure 8 shows the pressure
coefficient distribution and the velocity magnitude in lattice form around the cylinder defined at Re= 140000. Based on the results obtained by ELBM, flow field can be computed with reasonableaccuracy and stability at very high Reynolds number using ELBM. It should be noted that thebounce back boundary condition might be affected on the accuracy in some case [7]. In that case,proper boundary condition based on kinetic theory must be used [15].

Fig. 8 Pressure coefficient distribution and absolute velocity value at Re=140000

7. Conclusion
We have implemented Entropic Lattice Boltzmann method in parallel operation using CUDA
Fortran programming language on the NUT’s GPGPU system. Entropic LBM is used to perform
unconditionally stable simulation, while parallel computation is used to speed up simulations andimprove efficiency of code for the future large scale engineering simulations. The parallel code is
aimed to be used for large scale and long term simulation in engineering field in complicatedgeometry at high Reynolds number fluid flow. Numerical accuracy has been tested on doubly
periodic shear layer flow and compared with the solution of the standard LBM. Further, we solved
flow past a cylinder at Re=140000 to show performance of the ELBM parallel code. To evaluate
accuracy, Strouhal number at Re=1000, 3900 are defined and compared with experiment and othernumerical simulation results. The results were in good agreement. The parallel code of the ELBM is7

Transactions on GIGAKU 4(1)(2017)04006/1-8

implemented in the combination of low memory usage and fast computation. Using GPGPU for
ELBM, the computation can be performed 10 times faster than serial code. Based on the stable andaccurate computation of ELBM, parallel computation will be very efficient method for CFD.
Acknowledgments
Authors would like to express their gratitude to Prof. Tomohiro Degawa for whose valuableadvices and comments during the implementation and improvement stages of the parallel code.
References[1] Succi S. The lattice Boltzmann equation: for fluid dynamics and beyond. Oxford university

press, 2001.
[2] Bailey, Peter, et al. "Accelerating lattice Boltzmann fluid flow simulations using graphics
processors." 2009 International Conference on Parallel Processing.IEEE, 2009.
[3] Mittal, R., and S. Balachandar. "Effect of three‐dimensionality on the lift and drag of
nominally two‐dimensional cylinders." Physics of Fluids, Vol.7.8, pp. 1841-1865.1995[4] Karabelas, S. J. "Large eddy simulation of high-Reynolds number flow past a rotating
cylinder." International journal of heat and fluid flow, Vol. 31.4, pp. 518-527.2010.
[5] Rahman, MdMahbubar, MdMashud Karim, and Md Abdul Alim. "Numerical investigationof unsteady flow past a circular cylinder using 2-D finite volume method." Journal of
Naval Architecture and Marine Engineering, Vol.4.1, pp. 27-42.2007.
[6] Rajani, B. N., A. Kandasamy, and SekharMajumdar. "LES of Flow past Circular Cylinder at
Re= 3900." Journal of Applied Fluid Mechanics, Vol.9.3, 2016.
[7] Karlin, I. V., F. Bösch, and S. S. Chikatamarla. "Gibbs' principle for the lattice-kinetic
theory of fluid dynamics." Physical Review E, Vol.90.3, pp. 031302.2014.
[8] Ansumali, Santosh, et al. "Entropic lattice Boltzmann simulation of the flow past squarecylinder." International Journal of Modern Physics C, Vol.15.03, pp. 435-445. 2004.
[9] Ruetsch, G. and Massimiliano, F. "CUDA Fortran for scientists and engineers." NVIDIA
Corporation, 2701. 2011.
[10] Brown, David L. "Performance of under-resolved two-dimensional incompressible flow
simulations." Journal of Computational Physics, Vol 122.1, pp. 165-183.1995.
[11] Hou, Shuling, et al. "A lattice Boltzmann subgrid model for high Reynolds number
flows." arXiv preprint comp-gas/9401004. 1994.
[12] Liu, Tiancheng, et al. "Extended lattice Boltzmann equation for simulation of flows aroundbluff bodies in high Reynolds number." International Colloquim on Bluff Bodies
Aerodynamics, pp. 1-16. 2008.
[13] L. Ong, J. Wallace. The velocity field of the turbulent very near wake of a circular cylinder,
Exp. Fluids, Vol 20, pp. 441, 1996.
[14] Ayurzana, B. "Free surface flow through porous media by Lattice Boltzmann
method."Proceedings of the Mongolian Academy of Sciences, vol 56.01, pp. 42-54. 2016.
[15] Chikatamarla, S.S. and Karlin, I.V. "Entropic lattice Boltzmann method for turbulent flow
simulations: Boundary conditions." Physica A: Statistical Mechanics and its
Applications, 392(9), pp.1925-1930 , 2013
(Received: 7 October 2016, Accepted: 11 January 2017)8

Transactions on GIGAKU 4(1)(2017)04007/1-7

A Study on the Relationship between the Mindset of Taxi Drivers andthe Attribution of Accident Causes and Driving Behavior
Makiko Okamoto1, Katsuko Nakahira,T2,*
Department of System Safety, Nagaoka University of Technology,2)
Department of Information and Management System Engineering, Nagaoka University ofTechnology,1603-1 Kamitomioka-machi, Nagaoka 940-2188, Japan*E-mail: okamoto@kjs.nagaokaut.ac.jp
1)

In the field of Education, a learner’s mindset(the incremental mindset[abilities can be developed through effort] and the entity mindset [abilities are
static and unchangeable]) is thought to influence one’s learning motivation.
The prior enables learning motivation to be maintained even after failures
because the failure is attributed to a lack of effort, and the issue can thus beresolved with more effort. It has been shown that drivers who attribute theiraccidents to a lack of effort, experience fewer accidents. The purpose of thisstudy was to reveal the relationship between the mindset of drivers andattribution of accident cause, driving behavior. Results showed that there
were positive correlations between incremental mindset and the causal
attribution of “everyday effort”, higher accident rate and longer driving
distance. The group with incremental mindset(ITG) was more likely to haveaccident models that involved other persons. The study showed thepossibility that ITG had higher accident rates since they consideredthemselves to be able to “generate higher revenue,” rather than “prevent
accidents” with effort, and engaged in long-distance driving that could lead to
accidents. The results may be utilized in basic materials that can be provided
to drivers for precise safety education.

1. Introduction
It has been indicated that, in traffic safety education, it is desirable to teach accident prevention
in which the individual is the main factor in causing the accident, without attributing it to outside
factors such as others of traffic [1]. Using a questionnaire, Yamashita[2] found that, in a situation inwhich the occurrence of an accident is assumed, individuals who have an attitude of strictly taking
responsibility on themselves have a more desirable attitude toward safety, including egocentricity
and impulsivity, when compared to individuals who pass the responsibility onto others. Further, aquestionnaire survey revealed that the none accident group demonstrated higher attribution of thecause of the accident to themselves when compared to the accident-prone group. Yamashita[2]attempted to explain this as a causal attribution for “failure.” In the event of a failure, i.e., meetingwith an accident, the accident-prone group attributed it, not to themselves, but to the externalsituation or a mistake of the other party. In contrast, the low-accidents group attributed it to theirown lack of driving ability or lack of effort in confirmation, and as such, they have an increased
motivation to try to improve themselves, which makes them more safe drivers. Therefore, an attitude
of attributing responsibility to oneself is considered important in accident prevention.
On the other hand, the aforementioned research did not reveal which attribution to effort orability was important among the self-related factors. In this regard, it is meaningless whether it iseffort or ability, because, while the feeling that one is able to control the target of attribution isimportant, it has been pointed out that a distinction between the two is meaningless [3].
However, in terms of education, regarding the difference between causal attribution toeffort/ability, it is known that the student’s own ability combined with their perception of
intelligence (mindset), governs the pattern of behaviors and emotions that occur, mainly in the1Transactions on GIGAKU 4(1)(2017)04007/1-7

context of academic achievement. In this context, we need to consider the incremental theory andentity theory. The former is the idea that one’s ability changes, and if one puts forth effort, abilities
can improve. The incremental theory group (ITG), which focuses on the improvement of one’s own
ability and is aimed at the acquisition of knowledge (Mastery goal), attributes the cause of facing
difficulties and failings to personal lack of effort, is motivated by resolution using strategies different
from those of the past, and is less likely to experience negative feelings. On the other hand, the entity
theory is the idea that ability is innate and that it is difficult to control it. Therefore, the entity theory
group (ETG), which focuses on avoiding low assessments of one’s own abilities (performance goal),
attributes the cause of failure to personal lack of ability, and, as its adherents are of the position that
failure has exposed their own low capacity to those around them, negative feelings are likely tooccur, their motivation for subsequent challenges is reduced, and they become evasive [4]. Thus, inthe context of learning, the incremental theory is tied to an attitude of active participation in theclassroom and is generally considered more effective in the promotion of learning as comparedto the entity theory of intelligence [5]. Additionally, in the industrial sector, employees who
focus on the learning objective related to the incremental theory show a positive correlation
between effort and business performance, which is mediated by self-regulation relating tosetting goals, effort, and planning [6]. Therefore, a correlation between mindset andacademic/business performance is evident to some degree.
On the other hand, in the case of driving businesses like buses and taxis, not only
improvement of business performance, but preserving the safety of passengers and other people
in traffic is also important. However, the relationship between mindset and accidents that may
occur accompanying this business is not clear. Accordingly, the aim of this study was to reveal
the relationship between the driver’s mindset and his/her causal attributions of accidents, as wellas the driving behavior and accidents observed in professional drivers. The aim was to extract
points for developing an effective educational intervention.
The study was conducted with the ethics approval of Nagaoka University of Technology.

Survey 1
2-1 Experimental method(1) Mindset assessment scale (explicit mindset)
A questionnaire was designed based on a scale created by Hong et al. [7]. This questionnaire
contained 3 questions, and used six-point scale. In our survey, a higher response score indicates atendency to subscribe to the incremental theory. Note that, in order to clarify that “ability” pertained
to driving, “ability to drive a car” was used rather than “talent.” As this questionnaire attempted tomeasure the mindset that the experiment participants themselves were explicitly aware of, it wasreferred to as the explicit mindset.
(2) Causal attribution scaleWith reference to the work of Nasu et.al.[8], a self-assessment with the scene imagination
method was used. Specifically, subjects were asked to envision a scene in which they encountered an
accident while they were driving. The attributed cause of the accident had six classifications, theefforts just before (the efforts just before the accident directed towards avoiding it), normal effort(efforts directed at driving safely in normal circumstances), ability, physical condition, luck, andother traffic. Each item had two questions, with a total of twelve questions. Responses were made ona six-point scale. A higher value for a response demonstrates a higher amount of attribution to thatfactor.
(3) Mindset IAT (implicit mindset)
Mindset was measured using a self-report questionnaire that is used to assess an explicit mindset
obtained as a result of the individual’s introspection. However, as responses may be biased towards
social desirability and are influenced by personal introspection, it has been indicated that it may bedifficult to measure genuine mindset using this method. Accordingly, in addition to the questionnaire,
an implicit attitudes measure, known as the Implicit Association Test (IAT) [9], was used in thepresent study. The IAT is a method of measuring the strength of ties between concepts from the2

Transactions on GIGAKU 4(1)(2017)04007/1-7

response time taken while sorting words into specific categories in a vocabulary grouping task. It has
been demonstrated that unconscious behaviors conventionally considered difficult to measure can bemeasured to some degree by using this technique [10]. In the present study as well, we attempted toobserve correlations relating to accidents that could not be measured using an explicit indicator likea questionnaire.
IAT categories and stimulus words were created based on the work of Fujii et.al. [11]. Subjects
classified stimulus words associated with the attribution categories of “changes” and “does notchange,” and the two types of attribution categories of “intelligence” and “luck” in seven reaction
blocks into the presented categories, as fast as possible, using a key on a notebook PC. Note that,though the stimulus words were fundamentally derived from those used by Fujii at.al[11], as thepresent survey target was driving-related ability, “scholarly ability” was replaced with “driving
ability.” Scoring was conducted according to the technique suggested by Greenwald et al. [12]. A
higher score indicated a higher implicit tendency for subscribing to the incremental theory.

Intelligence

Luck

Change

Does not change

Driving ability

Figure 1 Example of IAT screen used in this survey
2-2

Experiment participantsParticipants were thirty-eight taxi drivers (all male, average age 57.34 years, SD=11.84 years) .
2-3 Procedure
The participants were led to an experiment hall prepared by the subcontracting company inadvance, and the details of the experiment were explained. It was conveyed that the experimentresults would be statistically analyzed, individual data would not be publicly announced
externally—including to the subcontracting company—in a format that could be identified, theexperiment results would not be related to professional assessment of the experiment participant, andthe experiment could be stopped midway if requested. Subsequently, after explaining the mindset
IAT implementation method while demonstrating it on the screen of the notebook PC, theexperiment was carried out, and when completed, the experimenter informed the participant
accordingly. Next, a questionnaire with the mindset assessment scale and the causal attribution scale
was presented in random order. After completing the questionnaire, the structure of the mindset IAT
and the research objectives were explained to the participant. Further, in order to confirm thetest-retest reliability of the mindset IAT, it was carried out a second time either one week or ten days
after the initial assessment. None of the participants reported that they wished to stop the experimentmidway.
2-4 Accident rate and mindset
We examined the relationship between the number of accidents occurring between December 16,
2009 and April 30, 2014, and mindset. Note that, as the period of membership with thesubcontracting company differed across participants, the number of accidents was divided by theinvestigation period to calculate the accident rate during the investigation period. For participants
who began membership with the subcontractor after the investigation period had begun, the numberof accidents was divided by the actual membership period.
2-5 Result
A positive correlation was observed between explicit mindset and normal effort (r=0.34, p<0.05).
A significant correlation between implicit mindset and causal attribution was not observed.
A significant correlation between explicit mindset and implicit mindset was not observed.
A relatively strong correlation was observed between the first and second mindset IAT results3

Transactions on GIGAKU 4(1)(2017)04007/1-7

Table1 Correlations of mindset, attribution of cause, accident rate

Explicit

Explicit
MS

Implicit
MS

Efforts
Just
BeforeNormalEffortNormalOthersOthersJust
Before

Mood,
PC

Luck

AbilityAccident
rate

1

0.2

-0.14

0.34*

-0.22

-0.19

-0.09

-0.21

0.11

0.50**

1

0.06

0.21

0.09

-0.1

-0.03

0.11

0.16

0.12

MS
ImplicitMS(r=0.48, p<0.05).
A positive correlation was observed between explicit mindset and accident rate (r=0.50, p<0.05),
whereas a significant correlation was not observed for implicit mindset(Table 1).

3.Survey 2
Distance traveled and mindset
The correlation between 7 months in 2014 distance traveled and mindset was investigated. Note
that, as this was an actual behavioral observation rather than a laboratory experiment; therefore, thesignificance level was set to 10%. Additionally, in order to control for the effect of the number of
workdays, the average distance traveled for one workday (the total kilometers traveled during theinvestigation period divided by the total word days) was used.
As a result, a weak positive correlation was observed between explicit mindset and average
distance traveled; a trend was observed in which the higher the mindset, the longer was the average
distance traveled (r=.28 p=.09). A significant correlation with implicit mindset was not observed.
3-2 Accident type and mindset
The experimenter extracted the participants’ behavior type at the time of an accident from theaccident circumstances and the sketch of the scene described in the accident reports collected by thesubcontracting company for the accidents that were included in the analysis for section 2-4, andcompared the findings between the incremental and entity groups. With regard to these groups, using
the average response on explicit mindset as a standard, those with a high average value were placedin the incremental group and those with a low value were placed in the entity group. As a result,when the incremental group was compared to the entity group, they had more accidents during
course changes, left and right turns, and departure. Meanwhile, the entity group had more
single-accidents while backing up and going straight (Table 2).
3-1

Table 2 Result of accident type and mindset
Behavior
Type
Group

Coursechanges

Left rightturns

Departure

Backingup

Goingstraight

Parking/s
topping

ITG

12

9

8

0

4

ETG

5

2

4

9

9

3-3

Total

11

Opening/
closing
doors
3

11

2

42

47

Accident report form entries and mindset
The accident circumstances were entered in the accident report form used in 3-1, by the driver
who was involved in the accident. This included the facts and causes of the accident. Consequently,
the cause of the accident was described based on their subjective judgment. Accordingly, from thedescription of the accident circumstances, the experimenter classified the participants based onwhether the cause of the accident was attributed to themselves or the other party. For example, if the4Transactions on GIGAKU 4(1)(2017)04007/1-7

participant used “...hit,” it would be attributed to themselves, and he used “...was hit by,” it would be
attributed to the other party. Accidents in which the causal attribution could not be understood byreading the description were omitted.
Using this data, when a chi-square test was conducted to compare the incremental and entity
groups in terms of their attributions of the cause of the accident, the incremental group was found to
attribute the cause to the other party more frequently (χ2=4.62, df=1, p<0.05)(Table 3).
Table 3 Result of causal attribution that was read from the description of the accident reports
and mindset
Group

ETG

ITG

Total

CausalThe Other Party

16

28

44

Attribution

Themselves

25

17

42

Total

41

45

86

4. Discussion
The results are summarized in Figure 2.
When those with ITG were asked to imagine a scene of encountering an accident and respond,
they generally attributed the cause of the accident to their own effort. However, when theyspecifically described the circumstances of an accident they had actually encountered, they exhibiteda tendency to consider the other party to be the cause. Additionally, in terms of behavior type, the
ITG had more accidents in settings where their own vehicle was accompanied with changes intraveling direction, like when changing course, turning left or right, or departing. In particular, whenlooking closely at the details of the descriptions of the circumstances in which the accidents occurred,
we observed more cases in which the course of the vehicle had changed according to the directionsof a passenger or when the driver found a passenger. Further, an overall tendency to travel for longer
distances for business was observed in the ITG.
As evident from the good performance of the employees with a strong learning goal focus that is
linked to the incremental theory, the concept that “ability can increase with effort” is thought to be
specifically represented in the idea that “if one puts forth effort, one’s performance improves.”
Considering that if one puts forth effort one’s business the revenue would increase, it is possible thatdriver believe that lengthening the distance traveled, actively looking for passengers, as well astrying to drive in response to the passenger’s wishes would help increase income.
Further, the tendency of ITG participants to attribute the cause to their own lack of effort when
imagining a general accident encounter, but to attribute the cause to other parties when they
imagine a scene of encountering an accident:attributed the cause to their own effort.

ITG

High accident
rate

they specifically described the circumstances of an
accident they had actually encountered,:attributed the cause to the other party
more accidents in settings where their own vehicle wasaccompanied with changes in traveling directionchanging course, turning left or right, or departing.

Tendency to travel for longer distances
for business

Figure 2 Summary of results

5Transactions on GIGAKU 4(1)(2017)04007/1-7

themselves are faced with an accident seems contradictory. However, taxi drivers have theopportunity to receive traffic safety education through lectures, etc. within their companies. At such
events, they are taught that their “efforts in preventing accidents are indispensable,” and as a result, it
is possible that they answered with the common belief that one’s own lack of effort is the cause of an
accident. Meanwhile, regarding cases in which they actually encounter an accident, the effect of aself-justification bias according to the cognitive dissonance theory needs to be considered [13].
Specifically, it is possible that when one encounters an accident while “one is putting forth effort” toimprove business revenue, he/she thinks “I, who was putting forth effort, caused an accident tooccur,” which leads to cognitive dissonance. Therefore, one may describe the other party as thecause to strike a balance.
Finally, we did not find a relationship between implicit mindset and driving behavior, etc.,
because all of the items surveyed were conscious behaviors, while the IAT used in this study
measured unconscious behavior [10]. As the descriptions on the accident report forms were filled inbased on the driver’s own explicit consciousness, i.e., distance was used as an indicator of
performance, and the driving behavior shown, such as turning left or right, was used as behavior
type at the time of the accident, all of which are consciously conducted behaviors, it is thought thatthey were not subjects for the IAT prediction.

5.Conclusion
The above findings indicate that the stronger the incremental tendency in the mindset of thedriver, the more they think, “if I put forth effort, my performance will increase,” and consequently,
their distance traveled increases, and their tendency to attach importance to the passenger becomesstronger. This may also lead to focusing too much on the presence of a passenger and not sufficiently
checking one’s surroundings, or driving continuously without breaks to increase performance.
Accordingly, with respect to the relationship between effort and performance, as an educational
intervention for this group of drivers, it is necessary to focus on helping drivers to be aware of theirown cognition, as well as having them understand that, though putting forth effort may increase
performance, more losses occur due to victim compensation, etc., if an accident is caused.
Additionally, they need to understand that it is therefore necessary to take appropriate breaks whiledriving and to consciously check one’s surroundings.
Additionally, as there is a tendency in ITG to̶while as a common belief considering the cause of
accidents to be one’s own lack of effort̶describe in the report that the cause of the accident wasanother party when actually faced with an accident, it is possible that there is a cognitive discrepancy
with respect to the common belief and the causal attribution at the time of an accident. When this is
not noticed by the persons themselves, the effectiveness of educational efforts aiming forunderstanding the accident as their own problem will likely also be insufficient. In terms ofeducation, prompting introspection about the cause of the accident one encounters and promoting
awareness of the aforementioned discrepancy are thought to be important.
Acknowledgments
This paper is obtained by revision and improvement on the basis of the manuscript of the next
conferences proceedings
T. Murohushi, K.T. Nakahira, M. Okamoto, “Study on assessment of participants’ theories ofintelligence: On the purpose of improvement in safety of car driving.” 12th Forum on
Information Technology, pp.545-548, 2013
M. Okamoto,”A study on the measure against improvement in safety consciousness of taxi drivers
based on safety culture” 28th International Congress of Applied Psychology, Paris France,
e-poster, 2014
References[1]国土交通省大臣官房運輸安全監理官室,”事故・ヒヤリハット情報の収集・活用の進め方〜
事故の再発防止・予防に向けて,” http://www.mlit.go.jp/common/001061869.pdf, April.
2009 (2016.10.3)
6

Transactions on GIGAKU 4(1)(2017)04007/1-7[2] N, A. Yamashita, “Study of the Measurement of the Safety Attitude,” IATSS Review, vol.16,
no.4, pp.17-24, 1990
[3] Y. Araki, “A review of causal attribution,” Socio-Environmental Studies, vol.5, pp.139-151, Mar.
2000.
[4] C.S Dweck,& E. L. Leggett, A social-cognitive approach to motivation and personality
Psychological Review 95, 256-273, 1988.
[5] C.S. Dweck, & A. Master, “Self-theories motivate self regulated learning, “ In: D.H. Schunk & B.
J. Zimmerman(Eds.) “Motivation and Self-regulated Learning : Theory, Research andApplications. NJ Lawrence Erlbaum Associates, pp.31-51, 2008.
[6] D.M. VandeWalle, S.P. Brown, W.L. Cron, J.W. Slocum, “The influence of goal orientation andself-regulation tactics on sales performance: A longitudinal field test.” Journal of Applied
psychology, no.84, pp.249-259, 1999
[7] Y, Hong, C. Chiu, C.S. Dweck, D.M. Lin, W.Wan, “Implicit theories, Attributions, and coping:
A meaning system approach.” Journal of Personality and Social Psychology, no.77,
pp.588-599, 1999
[8] M. Nasu, M. Horino, “Causal attributions and achievement-related affects.” Japanese Journal of
Educational Psychology, no.39, pp.332-340, 1991
[9] A.G. Greenwald, D.E,McGhee, J.L.K. Schwartz, “Measuring individual differences in implicit
cognition : The Implicit Association Test.” Journal of Personality and Social Psychology,
no.74, pp.1464-1480, 1998
[10] J.B. Asendorph, R. Banse, D. Mucke, “Duble dissociation between implicit and explicit
personality self-concept : The case of shy behaviour.” Journal of Personality and Social
Psychology, no.84,380-393[11] T. Fujii, H. Uebuchi,”Assessment of Participants’ Theories of Intelligence: Reliability andValidity of the Implicit Association Test.” Japanese Journal of Educational Psychology,
no.58, pp.263-274, 2010[12] Greenwald AG, Nosek BA, Banaji MR. Understanding and using the implicit association test :1. An improved scoring algorism. Journal of Personality and asocial Psychology, no.85, pp.
197-216, 2003.
[13] L. Festinger, “Informal social communication.” Psychological Review, no.57, pp.271-282, 1957
(Received: 7 October 2016, Accepted: 5 December 2016)7

Transactions on GIGAKU 4(1)(2017)04008/1-6

Engineering Design Exercise (EDE) for Faculty Development in
National Institute of Technology, Nagaoka College
Rumi Shindo Togashi 1,*, Taku Kiryu 2, Shigehiro Toyama 3, Yoshinori Tokoi 4,*, Tetsuro
Iyama 5, Fujio Ikeda 6, Yuki Murakami 7, Shin-ichi Akazawa 8, Yasuko Tsuchida 9
1)
Division of General Education, National Institute of Technology, Nagaoka College,
2)
Division of General Education, National Institute of Technology, Nagaoka College,3)
Department of Electronic Control Engineering, National Institute of Technology, Nagaoka
College,4)
Department of Electrical and Electronic Systems Engineering, National Institute of Technology,
Nagaoka College,5)
Department of Mechanical Engineering, National Institute of Technology, Nagaoka College,
6)
Department of Mechanical Engineering, National Institute of Technology, Nagaoka College,
7)
Department of Civil Engineering, National Institute of Technology, Nagaoka College,
8)
Department of Materials Engineering, National Institute of Technology, Nagaoka College,
9)
Division of General Education, National Institute of Technology, Nagaoka College,
888 Nishikatakai, Nagaoka, Niigata 940-8532, Japan*E-mail: rumi@nagaoka-ct.ac.jp

The National Institute of Technology, Nagaoka College (NITN)’s System
Design Innovation Center (SDIC) has launched several systems and tools
to encourage research promotion and education reform. Since 2015, SDIC
has applied engineering design exercise (EDE) to not only advanced
courses students but also new (or “rookie”) teachers, the latter of which
participate in those courses’ advanced experiments as part of their faculty
development (FD). This EDE focuses on four competencies: obtaining afundamental ability for cooperative working, understanding of innovativeeducational skills, building relationships among teachers from different
departments, and creating opportunities for regional contribution.
According to a distributed questionnaire to rookie teachers who
participated in the class as well as from gathered opinions of managers,
we can say that the understanding of active learning and relationshipdevelopment among new teachers did improve as a result of the EDE.
Some participants noted, however, that this program may require much
more effort than the organizers expected, so if this EDE is adopted byother institutes, the theme may need to be adjusted to suit the carrier aswell as the business volume of participants. We expect that this package
can be modified in terms of content to be suitable to each carrier, thereby
solving some potential problems.

1. Introduction
National Institute of Technology, Japan (KOSEN) has been required to improve and enhance its
education system and to improve the quality of its research. Various measures had previously beentaken; for example, the introduction of a model/core curriculum, a quality assurance program tooversee all education at KOSEN, among others. As one of these activities, faculty development (FD)for teachers and staff at the college has been introduced to encourage cooperation and coordination
between teachers and staff and to improve individual facilitation and instructional design skills. Eachcollege that is a part of KOSEN conducts lectures or trainings as part of FD, such as instruction
seminars at the school, study-abroad programs, and an exchange program between KOSENs. As
training for new teachers, some KOSENs have mentoring programs.
1

Transactions on GIGAKU 4(1)(2017)04008/1-6

The National Institute of Technology, Nagaoka College (NITN)’s System Design Innovation
Center (SDIC) has designed several systems and tools to encourage research promotion andeducational reform [1]. SDIC has also tried to enhance synergistic effects between active research andeducation. Since 2013, “advanced experiments” has been required for those students enrolled inadvanced courses. During the first semester of the year, engineering design exercise (EDE) wasarranged for students to learn how to manage the current rapid technological development andglobalization and to cope with recent changes in Japan’s industrial structure. EDE is designed toserve as a cross-field education program[2][3]. Within the EDE, instructors use teaching methods suchas active learning (AL) and project based learning (PBL).
Since 2015, SDIC has used EDE for not only advanced courses students but also for new (or“rookie”) teachers. It expects that the following four competencies will be achieved: (1) learn
fundamental competencies for collaborative working, (2) understand advanced examples of teaching
methods like AL and recognize essential skills or requirements for advanced knowledge and skills,(3) build relationships among teachers through cooperative research held both inside and outside of
the college, and (4) seek opportunities for regional contribution and community involvement. In thispaper, we propose and introduce EDE as a training method for rookie teachers. We also share thequestionnaire we developed for this group in 2015 and 2016, then discuss efforts, positive outcomes,
and problems that may remain to be solved.

Contents of Engineering Design Exercise (EDE) as Faculty Development
EDE is part of advanced courses students’ first semester as part of the required subject “advanced
experiments.” As one part of faculty development and induction, rookie teachers join and participate
in this class as students. This experimental attempt started in 2015. In the class, participants practicemethods needed for engineering design processes and solve real problems by using these methods.
Contents of the current curriculum—the purpose of which is mentioned in Section 1—are shown in
Table 1.
Table 1 Contents of Engineering Design Exercise

Contents
1. Guidance
Logical thinking- Ice-breaker session
2. Logical thinking- Meeting procedures
and consensus building

3. Logical thinking- Idea creation tools
and plan proposal

A. Need for engineering design exercise (EDE)
B. Ice-breaker session
A. Meeting procedures for engineering facilitation
B. Workshop (methods for creating new ideas)a. Brainstorming
b. KJ method
C. Methods for consensus buildinga. Votingb. Priority decision making matrix
c. Identify advantages and disadvantages
d. Payoff matrix
A. Idea creation (how to get to advanced and definite ideas)a. Logic tree (top-down and bottom-up)
b. MECE(mutually exclusive and collectively exhaustive)c. TRIZ (Teoriya Resheniya Izobretatelskikh Zadatch)d. Osborn's Checklist
B. Plan proposala. Process to create the plan
2

Classes
(150
min /
ea.)
1

1

1

Transactions on GIGAKU 4(1)(2017)04008/1-6

4. Logical thinking- Activities(plan proposal)

5. Engineering
facilitation- Presentation ofplanning
6. Engineering
facilitation andplanning for design
theme- Activities(improvement of plan)
7. Planning fordesign theme- Design review
8. Planning fordesign theme- Preparation forpresentation ofoutcomes
9. Presentation ofoutcomesb. Requirements for the outstanding plan
Practice methods for the meetings
Plan proposal about real problems
Preparation for the plan presentation<Examples (real problems)>
1. Development of a management system for compacting
concrete
2. Development of river flow velocity measurement techniques
3. Development of teaching materials using
LEGO MINDSTORMS® for technical courses in junior highschool(theme in 2015 and 2016 for rookie teachers)
Plan presentation and exchange of opinions
among students in the class, teachers and staff in the school,teachers in the junior high school, and supporting companies
Objective: Improvement of plan
Trial-based PDCA cycle (e.g., “plan, do, check, act”)
Improve the plan, based on opinions, comments, and potential
problems obtained at the plan presentation
Practical use PDCA cycle

Simple presentation and exchange of opinions between groups
dealing with the same theme
Decision on the final plan
Trial-based PDCA cycle
Prepare device, system, documents, etc. for result presentation

2

1

3

1

4

Results presentation
for students in the class, teachers and staff in the school, teachers
in the junior high school, and supporting companies
Evaluate the plan1

Investigation of the Effect of EDE as FDThis study used a questionnaire survey at the end of the in-class EDE in 2015 and 2016. The
questionnaire consisted of 13 questions that asked trainees about their feelings and perceived
difficulties throughout the EDE (Table 2).
Table 2 Questionnaire items
Q1. What was the proportion of work for EDE compared with your overall workflow?
Q2. Did you get a sense of the temperament of the Nagaoka KOSEN (NITN) students(i.e., advanced courses students) in EDE?
Q3. Do you think EDE should become one option for teachers’ faculty development?
Q4. The attempt to train new teachers along with students in EDE started in 2015. Do
you think it is a good idea to continue this next year?
Q5. Have the relationships among new teachers deepened through participation in
EDE?
3

Transactions on GIGAKU 4(1)(2017)04008/1-6

Q6. What do you think about the difficulty of this year’s design theme to develop
technology materials for a junior high school?
Q7. Did you have enough time in the class to carry out the design theme?
Q8. Did you need extra time to work on the design theme outside of your regular EDE
class?
Q9. Please tell me some good points and not-so-good points of EDE as FD.
Q10. In EDE class, the ways of logical thinking and engineering facilitation, such as
meeting procedures, consensus building, and plan proposal are presented. Have you
previously participated in a class or training session related to such facilitation skills?
Q11. In EDE class, what theme(s) do you see as being utilized in the future? (You canprovide multiple themes.)
Q12. Were you able to deepen your understanding of AL due to participation in the EDE
class?
Q13. EDE is supported by local companies, elementary schools, and junior high schools.
EDE education can provide measurable opportunities to contribute to the local
community. Was there any change in your opinion about the importance of this education
to contributing to the local community?
Q14. Please freely express your opinion or impression of EDE (e.g., did you notice
possible improvements to the theme that could be made next year?).
Questionnaire Survey Results and Discussion
The significant results of the questionnaire are shown in Figs 1 through 8. Fig. 1 reveals that all
trainees, both in 2015 and 2016, recognized the effect of EDE as FD, and in Fig. 2, trainees
answered that it was possible to interpret the students’ temperament. On the other hand, results in
Figs 3 through 8, which are concerned with the EDE program itself, are different when comparing
the years 2015 and 2016. The design theme in EDE class was the same in 2015 and 2016, so thisdiscrepancy probably came from the difference in the required effort of the majority of trainees. Fig.
7 shows that the effort needed for EDE varied among trainees and Fig. 9 shows that most trainees
were able to learn more about the AL method.
Moreover, the responses to Q9 (“Please tell me some good points and not-so-good points of EDE
as FD”) included the following. As the good points: “We can work with teachers who have different
professions,” “lectures on the conference method and consensus making were taught before tackling
the problems,” “we came to have further interest in students’ research and activities, knowing
first-hand what the students learn.” As for the not-so-good points: “It was absolutely impossible tofinish the activities for EDE within school hours.”
In addition, responses to the open-ended Q14 (“Please freely express your opinion or impressionof EDE”), one opinion was as follows: “It should be discussed whether the participation of rookie
teachers should be voluntary or mandatory.”
From the above, it is clear that participants were able to understand the temperaments of thestudents and understand the effects of AL as FD through the lectures about meeting method andconsensus making. While trainees considered that collaborating with teachers from different fields
was a key advantage, it was revealed that there was a difference between the trainees’ required effort.
In addition, only one of the trainees reported that the relationship among trainees worsened as the
EDE progressed. This may be due to the difference in trainees’ awareness regarding the importanceof their participation in EDE. Solving this problem will require further improvements to activity
efficiency.
As the background instruction of this EDE as FD is carried out, training or induction into a class4

Transactions on GIGAKU 4(1)(2017)04008/1-6

that already exists could reduce or, in some cases, eliminate the cost to invite more instructors.
Attending the lectures for students, however, in addition to their usual work may feel like excessive
effort to some trainees. Also, newly appointed teachers who are considered “general educators” maylack even a base level of engineering knowledge, which may result in more time required tounderstand the contents of EDE.

Fig. 1 Survey result (Q3)

Fig. 2 Survey result (Q2)

Fig. 3 Survey result (Q1)

Fig. 4 Survey result (Q4)

Fig. 5 Survey result (Q5)

Fig. 6 Survey result (Q6)

Fig. 7 Survey result (Q7)

Fig. 8 Survey result (Q8)5Transactions on GIGAKU 4(1)(2017)04008/1-6

Fig. 9 Survey result (Q12)Conclusion
EDE utilizes the required subject of advanced experiments as an FD opportunity for rookie
teachers. This attempt started in 2015 and aims to build fundamental competencies for collaborative
working, understanding and obtaining the skills required for those in higher education, building
relationships among teachers across different departments, and creating opportunities for regional
contribution. In an EDE class, logical thinking and engineering facilitation, such as meeting
procedures, consensus building, and plan proposal are presented. EDE participants try to solve real
problems using their skills and materials and methods presented in the lectures and then present theoutcomes. According to the results of the questionnaire provided to the participant rookie teachers, it
is clear that their understanding of active learning and the relationships among the new teachers
became better, while some considered that this program required much more effort than initially
expected. As such, if this EDE is adopted by other institutes, the theme should be customized to suitthe carrier as well as the business volume of participants. We expect that this package and system
can be improved and result in positive effects by modifying or customizing the content to be moresuitable to each carrier, solving their unique potential problems.
Acknowledgments
The authors wish to express gratitude to the supporting companies, participants andco-operators of the education component, rookie teachers who joined the class, colleagues in SDIC,
and NITN, the latter of which were crucial to permitting these trials in the first place.
References[1] Toyama, S., Ikeda, F., Iyama, T., Tokoi, Y., Akazawa, S., Murakami, Y. and Tsuchida, Y.,
System Design Education Program to Produce Innovative Personnel, Transactions of ISATE 2015,(2015) 277-280.
[2] Iyama, T., Toyama, S., and Tokoi, Y., A Practice of Engineering Design Education through theDevelopment of Educational Materials using LEGO-Mindstorms, The 4th International GIGAKU
Conference in Nagaoka, EO-17, (2015).
[3] Toyama, S., Miyakoshi, K., and Tasaki, Y., Design and Practice of Engineering Design
Education, Journal of Education in Japanese College of Technology, No.38, (2015) 67-72.
[4] Maeno, T., Official Guide Book to KEIO SDM: System and Design Thinking for Creating
Innovation. Tokyo: Nikkei BP, (2014).
[5] Oishi, K., Engineering Facilitation. Tokyo: Morikita shuppan, (2011).
[6] Pink, D.H., High Concept, the times of innovative personnel, translated by Kenichi, O., Tokyo:
Mikasa shobo, (2006).
(Received: 6 October 2016, Accepted: 6 December 2016)6