JSST2022

The 41st JSST Annual International Conferenceon Simulation Technology

Student Session Proceedings
Kyushu Institute of Technology, Iizuka, Fukuoka, JAPAN
Online Conference, Aug. 31—Sep. 2, 2022

Supporting CompaniesiiCooperating Societies

The Institute of Electrical Engineers of JapanThe Institute of Electrical Engineers of Japan, Kyushu Branch
The Japan Society for Industrial and Applied Mathematics
The Japan Society of Plasma Science and Nuclear Fusion Research
The Japan Society of Mechanical Engineers
The Visualization Society of JapanThe Robotics Society of Japan

iii

Session ScheduleDay 1: Wednesday, 31st August, 2022
Time

Hall A

Hall B

Hall CHall DHall E

Opening Ceremony (Hall A)

09:45 – 10:00
10:00 – 10:30

Plenary Talk 1 (Hall A)
10:30 – 11:00
11:00 – 11:15

Group Photo

11:15 – 11:30
11:30 – 11:45
11:45 – 12:00

Invited Talk 1
Tutorial Talk 1
Session 1 (2)OS 6 1
Session 6 (4)OS 7 1
Session 9 (3)12:00 – 12:15
12:15 – 12:30

Lunch Time

12:30 – 13:30
13:30 – 13:45
13:45 – 14:00

Corporate
Seminar 1A

Corporate
Seminar 1B

14:00 – 14:15
14:15 – 14:30
14:30 – 15:00

Plenary Talk 3 (Hall A)
15:00 – 15:30
15:30 – 15:45
15:45 – 16:00
16:00 – 16:15
16:15 – 16:30

Symposium 1 1
Session 2 (3)OS 8 1
Session 4 (3)OS 2 1
Session 7 (3)Symposium 1 2
Session 3 (2)OS 8 2
Session 5 (2)OS 2 2
Session 8 (2)OS 7 2
Session 10 (3)OS 4 1
Session 12 (3)OS 7 3Session 11 (3)OS 4 2
Session 13 (3)16:30 – 16:45
16:45 – 17:00
17:00 – 17:15
17:15 – 17:30
17:30 – 17:45
17:45 – 18:00

iv

Day 2: Thursday, 1st September, 2022
Time

Hall A

Hall B

Hall CHall DHall E

Symposium 1 3
Session 14 (3)Symposium 2 1
Session 17 (3)Invited Talk 3
Tutorial Talk 3
Session 20 (2)OS 7 4
Session 23 (3)OS 4 3Session 26 (3)Symposium 3 1
Session 21 (3)OS 5 1
Session 24 (3)Symposium 3 2
Session 22 (2)OS 5 2
Session 25 (2)09:00 – 09:15
09:15 – 09:30
09:30 – 09:45
09:45 – 10:00
10:00 – 10:15
10:15 – 10:30
10:30 – 10:45
10:45 – 11:00

Symposium 1 4
Session 15 (3)Symposium 2 2
Session 18 (2)OS 4 4
Session 27 (2)11:00 – 11:15
11:15 – 11:30
11:30 – 11:45
11:45 – 12:00

Symposium 1 5
Session 16 (2)12:00 – 12:15

Symposium 2 3
Session 19 (3)12:15 – 12:30

Lunch Time

12:30 – 13:30
13:30 – 13:45
13:45 – 14:00

Corporate
Seminar 2A

Corporate
Seminar 2B

14:00 – 14:15
14:15 – 14:30
14:30 – 14:35

Announcement about JSST2023 (Hall C)

14:35 – 15:15

Shotgun Session (Hall C)

15:15 – 15:30
15:30 – 16:30

Student Session + Corporate Exhibition (Remo)
16:30 – 18:00

v

OS 7 5
Session 28 (3)Day 3: Friday, 2nd September, 2022
Time

Hall A

Hall B

Hall CHall DHall EOS 1 1
Session 38 (3)OS 9 + OS 10 1
Session 41 (3)OS 1 2
Session 39 (3)OS 9 + OS 10 2
Session 42 (3)09:00 – 09:30

Plenary Talk 2 (Hall B)
09:30 – 10:00
10:00 – 10:15
10:15 – 10:30
10:30 – 10:45
10:45 – 11:00

OS 3 1
Session 29 (3)Invited Talk 2
Tutorial Talk 2
Session 32 (2)Symposium 3 3
Session 35 (3)Symposium 2 4
Session 33 (2)Symposium 3 4
Session 36 (2)11:00 – 11:15
11:15 – 11:30
11:30 – 11:45
11:45 – 12:00
12:00 – 12:15

OS 3 2
Session 30 (3)12:15 – 12:30

Lunch Time

12:30 – 13:30
13:30 – 13:45
13:45 – 14:00

OS 3 3Session 31 (2)Symposium 2 5
Session 34 (2)Symposium 3 5
Session 37 (2)OS 1 3Session 40 (2)

14:00 – 14:15
14:15 – 14:30
14:30 – 15:00

Closing Ceremony (Hall A)viSymposium1
Plenary Talk
Presenter:

Prof. Toshihiro Hanawa
Information Technology Center, The University of Toky

Title:
Smart Supercomputing by Integration of Simulation, Data, and Learnin

Abstract:
Traditionally, supercomputers have been used for Computational Science andEngineering (CSE) simulations. For future supercomputing, we should
consider integrating data analysis and machine learning to leverageknowledge in a smarter way than just traditional simulation.
This talk will introduce our activities toward “smart supercomputing,” including
Wisteria/BDEC-01 in the Information Technology Center, the University ofTokyo, and the “mdx” data platform under the joint operation with nine
national universities and two national research institutes in Japan

Biography:
Toshihiro Hanawa received his M.E. degree and his Ph.D. degree incomputer science from Keio University in 1995 and 1998. He was an
assistant professor at Tokyo University of Technology, Japan, from 1998 to2007, a research fellow at Center for Computational Sciences (CCS),
University of Tsukuba, from 2007 to 2008, and an associate professor of CCS
from 2008 to 2013. He was a project associate professor at Information
Technology Center, The University of Tokyo, from 2013 to 2015 and an
associate professor from 2015 to 2020. Since Dec. 2020, he has been aprofessor at the Information Technology Center, The University of Tokyo.
He is also the division leader of the Operations Support Division in JCAHPC(Joint Center for Advanced High-Performance Computing), a joint
organization by the University of Tsukuba and the University of Tokyo tooperate the Oakforest-PACS supercomputer system. He leads the design ofvii

the next Oakforest-PACS II system.
His research interests include computer architecture, interconnection
networks, and acceleration using GPU and FPGA

viii

Invited Talk
Presenter:

Prof. Mitsuteru Asai
Department of Civil and Structural Engineering,
Faculty of Engineering, Kyushu Universit

Title:
Multi-scale and -physics particle simulations for natural hazar

Abstract:
Particle methods such as the SPH and MPS methods can compute more
ef ciently problems involving large deformation, splitting, and coupling of theobject of interest than mesh-based methods such as the nite element
methods. For example, when applied to the design of disaster prevention andcountermeasures, free-surface ow problems in complex geometries, such as
tsunami runup and ooding in heavy rain disaster, can be easily solved, andeven multi-scale and -physics phenomena involving the failure of solids, suchas seepage collapse and scouring of breakwaters and seawalls. This talk will
share the recent progress on the multi-physics particle simulation related tonatural disaster issues

Biography:
Mitsuteru Asai received the M.E. and Dr. Eng. degrees from Tohoku
University in 2000 and 2003. After working at the Ohio State University as a
postdoctoral fellow, he became an Asistant Professor at Ritsumeikan
University in 2005. Since 2007, he has been an Associate Professor in the
Department of Civil Engineering, Kyushu University. His research area iscomputational mechanics, computational science, and VR/AR related tonatural disaster issuesixTutorial Talk
Presenter:

Dr. Shintaro Kawahara
Research Institute for Value-Added-Information
Generation (VAiG), Japan Agency for Marine-Earth
Science and Technology (JAMSTEC

Title:
VisAssets: A Visualization Framework for Unit

Abstract:
This tutorial introduces the procedure for developing visualization applications
using VisAssets, a framework that runs on the Unity game engine. VisAssets
provides the elements that constitute a general visualization ow as module
icons and visualization applications with GUI can be built by simply
connecting them appropriately on Unity Editor. This tutorial will show you
examples of building applications to visualize scalar and vector eld datausing VisAssets
VisAssets can download from the following URL.
https://github.com/kawaharas/VisAssets

Biography:
Shintaro Kawahara received Dr. Eng. from Ibaraki University in 2001. After
working as a exible researcher at the Ibaraki Prefectural Industrial
Technology Center and as an assistant at the Information Processing Centerof Ibaraki University, he joined the Japan Agency for Marine-Earth Scienceand Technology (JAMSTEC) in 2003. He belonged to the Earth Simulator
Center and was engaged in research on visualization of large-scale
simulation data and observation data. He is also involved in the disseminationof visualization software that he has developedxSymposium2
Plenary Talk
Presenter:

Prof. Kensei Kobayashi
Department of Chemistry Presentation,
Yokohama National Universit

Title:
Origins of Bioorganic Compounds and Their Homochirality: Approaches from
Simulation Experiment

Abstract:
How life was generated on the Earth? Are there any extraterrestrial bodies
that foster life on them? The novel scienti c eld to challenge such questions
is astrobiology. Laboratory simulation is one of the most useful approaches
to study prebiotic chemical evolution. In this talk, I will present
simulation experiments to synthesize bioorganic molecules like amino acids
in early Earth or extraterrestrial environments. It was proved that high-energy
particles irradiation was quite effective to produce them. Amino acids arechiral molecules and we utilize mostly L-isomers in our biosystem.
Simulation experiments to prove the origin of the homochirality of
biomolecules will also be introduced, where such quantum beams as
circularly polarized light and spin-polarized leptons were used

Biography:
Kensei Kobayashi is currently a professor emeritus of Yokohama
National University. He received his D. Sc. In Department of Chemistry,
University of Tokyo in 1982. In 1982-86, he studied in Laboratory of Chemical
Evolution, University of Maryland, U.S.A. He has been working in Yokohama
National University (YNU) since 1987, and was a professor at YNU in
2003-2020. He has been studying astrobiology, especially focused on
prebiotic chemistry toward generation of life and detection of life in extreme
terrestrial and extraterrestrial environmentsxi

Invited Talk
Presenter:

Mr. Ryuta Kawanami
Graduate School of Science and Technology,
Kyoto Institute of Technolog

Title:
Molecular dynamics study of microscopic mechanism of
OH radical-induced DNA damag

Abstract:
DNA sugar-phosphate backbone has sugars joined by phosphodiester bonds.
Highly reactive hydroxyl radicals are produced by ionizing radiation in a cell,which can abstract any hydrogen atoms from DNA. Carbon centered sugar
radicals produced by the hydrogen abstraction triggers a common lesion,
single-strand break (SSB) which induces disease such as carcinogenesis.
However, the abstraction mechanisms at the molecular level remain
unsolved. Thus, we employed molecular dynamics (MD) simulations onhydroxyl radicals around DNA to evaluate their spatial distribution andaccessibility to DNA. In particular, we introduced restraining potentials on thehydroxyl radicals, which enabled us to suf ciently sample the hydroxyl
radicals around the sugar moiety of DNA. The result showed that theaccessibility of hydroxyl radicals to each hydrogen atom on DNA was very
different from the evaluations reported previously. This indicates that theprevious evaluations based on the solvent-accessible surface area, or thewater accessibility are needed to be modi ed

Biography:
Ryuta Kawanami is a doctoral student in soft materials chemistry at Kyoto
Institute of Technology. He is interested in physical mechanisms of soft matterin lifexii

Invited Talk
Presenter:

Prof. Naohisa Yahagi
Graduate School of Media and Governance,
Keio Universit

Title:The Secret to Successful medical-industrial Collaboratio

Abstract:
The eld of medicine is often described as a complex system. However, it is
not dif cult for physicians to make diagnostic and therapeutic decisions basedon their perception of a patient’s condition and its changes to nd the data
representing the pathology accurately
In other words, the key to successful medical-engineering collaboration
begins with the engineers accurately understanding the tacit knowledge of thephysician who knows the patient well, while the physicians accurately
understand the advanced technology and theory. This interrelationshipenables the development of high-precision technology
In this lecture, I would like to talk about a case in which such an overly
complicated event was derived from simple medical common sense,
dreaming of the day when the solid results of medical engineeringcollaboration with you will surpass the world

Biography:
Born in Palo Alto, the USA, in 1974. I graduated from Keio University Schoolof Medicine and completed the doctoral program at Keio University Graduate
School. Pediatrician. Completed (summa cum laude) in the medical businessand management course at the University of Tokyo. Invited to and completed
the Harvard Business School “Managing healthcare delivery” as a
Scholarship. Designed and launched the world’s rst Clinical DataManagement Network, enabling centralized management of medical
information throughout Japan. Engaged in research and development of
automated diagnosis and treatment support technology based on diseasestate change prediction. Head of Digital Health, Digital Agency
xiii

Symposium
Plenary Talk
Presenter:

Prof. Shinobu Yoshimura
Department of Systems Innovation,
School of Engineering, The University of Toky

Title:
Super-simulations of Offshore Wind Farm on Fugak

Abstract:
An offshore wind farm consists of tens of large-scale wind turbines. To
promote such offshore wind farms in Japan, we need to consider severe
weather conditions caused due to typhoons and narrower sites. To do so, weare developing multiscale and multiphysics integrated simulations for an
entire offshore wind farm. The simulations consist of the following four
components, (a) Large-scale LES simulator of a whole wind farm considering
weather and terrain, RC HPC, (b) Large-scale LES simulator of two tandem
wind turbines to evaluate an effect of wake into power generation, FFR, (c)ow-induced vibration of wind turbine, ADVENTURE_Solid and
ADVENTURE_Fatigue, and (d) a parallel coupling tool to deal with FluidStructure Interaction, REVOCAP_Coupler. In this presentation, we describe
Fugaku-based super-simulations of an actual German offshore wind warm in
the North see, Alpha Ventus, consisting of twelve 5MW wind turbines

Biography:
Prof. S. Yoshimura is currently Vice President and Professor, The Universityof Tokyo (UTokyo). In 1987, he completed his studies in UTokyo (specialized
in nuclear engineering) as Doctor of Engineering. He became Lecturer,
Associate Professor in UTokyo, then Professor since 1999. His specialties are
High-performance and Intelligent Computational Mechanics with Real World’s
Applications. He serves as a Vice President of IACM and the President of
APACM, and has received numerous awards, including IACM Fellow Award(2014) and APACM Computational Mechanics Award (2013)xiv

Invited Talk
Presenter:

Dr. Prakasha Chigahalli Ramegowda
Institute of Microelectronics (IME)/
Sensors, Actuators & Microsystems (SAM),
Agency for Science, Technology and Research(A*STAR), Singapor

Title:
Coupled multiphysics simulation model for analyzing piezoelectric sensors
and actuator

Abstract:
The key to the design of advanced micro electro–mechanical (MEMS)
piezoelectric devices for actuator and sensor applications surrounded by uid
media is an accurate evaluation of the electromechanical, uid and circuit
interaction leads to a general formulation as a uid-structure–piezoelectric–
circuit interaction. These four- eld coupled problems can be decomposed into
the uid-structure, structure–piezoelectric, and piezoelectric–circuit
interactions using hierarchical decomposition and exchange of the variables
using partitioned approach. In this talk, we present computational methods[1,2,3,4,5] to solve structure–piezoelectric interaction, piezoelectric–circuit
interaction, and uid-structure-piezoelectric interactions in piezoelectricbimorph actuators and energy harvester
References[1] P. C. Ramegowda, D. Ishihara et al., International Journal ofComputational Methods, Vol.16, 1850106, 2019.
[2] P. C. Ramegowda, D. Ishihara et al., Finite Elements in Analysis andDesign, Vol.159, pp. 33-49, 2019.
[3] P. C. Ramegowda, D. Ishihara et al., Computer Methods in Applied
Mechanics and Engineering, Vol. 365, 113006, 2020.
[4] P. C. Ramegowda, D. Ishihara et al., Composite Structures, Vol. 245,
112284, 2020.
[5] D. Ishihara et al., Computers & Structures, Vol. 253, 106572, 202

xv

Biography:
Dr. Prakasha Chigahalli Ramegowda received a BE in Mechanicalengineering from the National Institute of Engineering in Mysore, India, in2010. He was with the MEMS design group at Society for Integrated Circuit
technology and Applied Research, India, as an Engineer, from 2011- 2015. In
2019 he completed his PhD in Mechanical Engineering from Kyushu Instituteof technology, Japan. He is a recipient of the prestigious JSPS Postdoctoral
Fellowship (Standard FY2019-22) from the Japan Society for Promotion of
Science. He is currently Scientist at Institute of Microelectronics (IME),
Agency for Science, Technology and Research (A*STAR), SingaporexviTutorial Talk
Presenter:

Dr. Hiroshi Watanabe
DEP Systems Business Unit,
Technical Product & Support Promotion, Hexago

Title:
Fusion of CAE and 3D Metrolog

Abstract:
Recently, we often hear the term digital twin, but it is used in a wide range
concept. It can be said that the common term is to reconstruct the event in
the cyber physical space based on the digital data obtained by observing theevent in the real world, and to predict the future by simulation. Some people
claim that the traditional CAE is a digital twin. This is based on the fact that a
product is built based on CAD, which is a blueprint, and the product ismanufactured according to it. However, as is well known, it is dif cult toprocess a product according to a blueprint. The accuracy of the simulation
can be said “dizygotic twins” (similar because they are born from the same
parent). In this paper, we will consider the points for constructing a digital twin
suitable for “identical twins” in the eld of computational engineering

Biography:
Manufacturing Intelligence division, Hexago
Vice President, JSCES

xvii

Committee
General Chair:
Hiroaki Nakamura (National Institute for Fusion Science)
Susumu Fujiwara (Kyoto Institute of Technology
Conference Chair:
Daisuke Ishihara (Kyushu Institute of Technology
Conference Co-Chair:
Katsuhiko Yamaguchi (Fukushima University)
Amane Takei (University of Miyazaki
Program Chair:
Seiki Saito (Yamagata University)
Hiroshi Tamura (Chuo University
Program Co-Chair:
Ryosuke Ueda (Tohoku University)
Yoshihide Shibata (National Institute of Technology, Gifu College)
Yoshihisa Fujita (Ritsumeikan University)
Kazuhiro Fujita (Saitama Institute of Technology
Post-Conference Chair:
Hiroaki Ohtani (National Institute for Fusion Science
Post-Conference Co-Chair:
Tomoko Mizuguchi (Kyoto Institute of Technology)
Susumu Fujiwara (Kyoto Institute of Technology)
Keisuke Nakano (Niigata University
Publication Chair:
Taku Itoh (Nihon University)

Publication Co-Chair:
Soichiro Ikuno (Tokyo University of Technology)
Hiroto Tadano (University of Tsukuba)

xviii

Finance Chair:
Tomoya Niho (Kyushu Institute of Technology
Publicity Chair:
Nobuaki Ohno (University of Hyogo
Exhibition Chair:
Amane Takei (University of Miyazaki
Conference Adviser:
Takuo Yasunaga (Kyushu Institute of Technology)
Takahiro Kenmotsu (Doshisha University
Committee Member:
Takashi Hara (National Institute of Technology, Kagoshima College)
Masafumi Yoshida (National Institute of Technology, Ube College)
Yuki Uchida (National Institute of Technology, Nagaoka College)
Makoto Nakamura (National Institute of Technology, Kushiro College

xixTable of Contents for Student Session
Day 2: Thursday, 1st September, 2022
Student Session (Poster Session)
Time: 15:30 – 18:00
Room: Remo
Chair: Hiroshi Tamura (Chuo University) andYoshihide Shibata (National Institute of Technology, Gifu college)P1-01 Explosive crystallization of amorphous Ge studied by molecular-dynamics
simulations ..........................................................................................1
Shunya Nagaoka and Manabu Ishimaru
P1-02 Effect of voltage application rate on resistance change characteristics in Ta2 O5 -based
ReRAM cells ....................................................................................... 5
Gaku Ishizu and Yusuke Nishi
P1-03 Simulation of Panic Buying in a Disaster Focusing on Consumer Information
Transfer ............................................................................................. 9
Taku Ideguchi, Kazuo Furuta and Taro Kanno
P1-04 Development of Multiscale Synchronic Simulation Method for Brain NeuronalGroups ............................................................................................. 13
Sho Ikeda
P1-05 Numerical prediction of water–ice phase–change with natural convection including density inversion region ............................................................................. 17
Kazuki Nishimoto, Daisuke Toriu and Satoru Ushijima
P1-06 Consideration of Dust Accumulation Process of Grows Fractal-wise ....................21
Masaki Ohtani and Norikane Kanai
P1-07 Base station placement problem in 5th generation mobile communication systems ... 25
Shoya Kurebayashi and Hiroshi Tamura
P1-08 Algorithm for Detour Hinge Vertex Problem of Circular-arc Graphs .................... 29
Tomonari Izumi, Syoma Nameki, Yoko Nakajima, Hirotoshi Honma and Toshihiro Fujito
P1-09 Development of a System to Predict Child Abuse in Newborns .......................... 33
Yuya Sato, Yoko Nakajima, Hirotoshi Honma, Yasuko Makita and Kazuyo Matsuura
P1-10 A System to Generate Kanji Reading and Writing Questions for Learners of the Japanese
Language .......................................................................................... 37
Toshiki Aisaka, Yoko Nakajima, Hirotoshi Honma and Tomoko WatanabexxP1-11 Molecular Dynamics Simulation of Clustered DNA Damage Composed of
Apurinic/Apyrimidinic Sites ....................................................................41
Kazushi Terakawa, Susumu Fujiwara, Tomoko Mizuguchi, Yoshiteru Yonetani, Naoya
Shikazono, Ken Akamatsu and Hiroaki Nakamura
P1-13 Automatic Chord Estimation Using CGAN ..................................................45
Ryuto Morita and Takashi Hara
P1-14 Identiﬁcation of Road Signs Using Shape Estimation from Depth Images ..............49
Haruka Aiko, Kazuhiro Takeda and Takashi Hara
P2-01 Enhancement of Images Reﬂected on Electronic Mirrors ..................................53
Shoma Higashiyama, Shota Furukawa, Takashi Hara and Kazuhiro Takeda
P2-02 Numerical Experiments of Swelling Objects Interacting with Newtonian Fluids ...... 57
Niku Guinea, Daisuke Toriu and Satoru Ushijima
P2-03 Improvement of Lagrange ﬂuid simulation for application to a linear plasma
experiment .........................................................................................61
Morihiro Mine and Kazuo HoshinoP2-04 Consideration of “Preventing deforestation and Forest preservation activities in
Myanmar” ......................................................................................... 65
Myo Thiri Thu and Akira Otsuka
P2-05 Analysis of Hydrogen Bonds inside Reverse Micelles by Molecular Dynamics
Simulations ........................................................................................70
Yuri Ono, Tomoko Mizuguchi and Mao Fukuyama
P2-06 Simulation of radiation source location estimation using Compton camera .............74
Yoshiki Ohori and Katsuhiko Yamaguchi
P2-07 Veriﬁcation of Estimation Accuracy for Shield Distribution Aggregation using Machine
Learning ........................................................................................... 78
Yuto Kondo, Masaharu Matsumoto, Kenji Suzuki and Katsuhiko Yamaguchi
P2-08 Electrostatic ﬁeld simulation using deep learning ...........................................82
Tsubasa Sagara, Masaharu Matsumoto, Kenji Suzuki and Katsuhiko Yamaguchi
P2-09 Multiscale analysis of magnetic materials in the mesoscopic region .....................86
Yasuyuki Igarashi and Katsuhiko Yamaguchi
P2-11 The Beneﬁt of Data Augmentation Performed by Surrogate Modeling ..................89
Takumi Nomura and Tomonori Yamada
P2-12 Conﬁgurational Entropy of a Monatomic Glassy System ..................................93
Ayata Ueno and Tomoko Mizuguchi
P2-13 Numerical analysis of hydrodynamic equations using Spreadsheet software ........... 96
Sakuya Inoue, Kiyomi Yamane, Munkhbileg Batchuluun and Hiroki Hommaxxi

P2-14 Analysis of cycle-to-cycle variation components in engine cylinder ﬂow ............. 100
Shizuka Morikawa and Yuichi Itohxxii

JSST2022

Explosive crystallization of amorphous Ge
studied by molecular-dynamics simulations
Shunya Nagaoka, Manabu Ishimaru
Department of Materials Science and Engineering, Kyushu Institute of Technology
nagaoka.shunya177@mail.kyutech.jp

Abstract. Anomalous rapid crystallization of several meter/sec, so-called explosive crystallization, often occurs during pulse laser irradiation and flash lamp annealing of amorphous IVgroup elemental semiconductors. To clarify the origin of rapid crystallization, crystallization
processes of amorphous Ge was examined by molecular-dynamics simulations using the Tersoff interatomic potential in combination with the Langevin thermostats. Amorphous Ge networks containing 125,000 atoms were annealed under a temperature gradient. It was found thatthe crystallization rapidly progresses from the high temperature to low temperature side, andlarge crystal grains were obtained. From the mean-square atomic displacement, it was suggested that a liquid layer plays an important role for inducing the explosive crystallization.

Keywords: Molecular-dynamics simulation, Germanium, Crystallization, Amorphous

1. Introduction
Crystallization of amorphous Si and Ge is one of the possible ways to fabricate polycrystalline
thin film transistors (TFTs). Several researchers have reported that abnormal rapid crystallization, so-called “explosive crystallization”, occurs by pulse-laser irradiation and flash-lamp annealing. The explosive crystallization was observed to proceed over 10 m/s [1]. It is thoughtthat the liquid layer at an amorphous/crystalline interface plays an important role for enhancing
the crystallization velocity. Since the explosive crystallization proceeds in a very short time, it
is difficult to experimentally clarity the origin of the rapid crystallization. Molecular-dynamics(MD) simulations can track atomic behavior with high spatial and temporal resolution, andtherefore they are useful for studying the crystallization processes.
We have previously studied the crystallization processes of amorphous Si in contact withcrystal Si during isothermal heat treatment by MD simulations [2]. On the other hand, MD
simulations of explosive crystallization require heat treatments of amorphous without crystallites under a temperature gradient. In the present study, we examined the crystallization processes of amorphous Ge under annealing with a temperature gradient.

1

JSST2022

2. Experimental
The crystallization processes of amorphous Ge were examined by MD simulations based on
the Tersoff3 interatomic potentials [3]. A cubic cell consisting of 25×25×25 diamond unit
cells contained 125,000 Ge atoms. The amorphous Ge was prepared by the melt-quench
method. The crystal Ge was first heated at 3500 K to produce liquid Ge. Then, liquid Ge wascooled to 300 K with a quenching speed of 1012 K/s, to obtain amorphous Ge. The amorphous
Ge obtained here was annealed at 2000-1700 K with a temperature gradient to observe crystallization from amorphous Ge. The temperature was controlled by Langevin method. The simulation was performed under a constant pressure condition and time step of 2.0 fs. MD simulations were carried out with an MD open-source code, LAMMPS [4]. The crystalline and amorphous regions were judged by using an Open Visualization Tool, OVITO [5,6].

3. Results and Discussion
Fig. 1(a) shows the structural evolution as a function of annealing time. To make it easy to see
the atomic arrangement, only the crystallized area is indicated. The blue and orange atoms
denote the cubic and hexagonal diamond structures, respectively. A temperature gradient was
set to be 1800 K at the top and 1000 K at the bottom. At 14.0 ns, a few crystal nuclei with adiamond structure preferentially appear near the high-temperature area. From 14.0 to 16.0 ns,
the crystal growth rapidly proceeds from high-temperature to low-temperature area. This rapid
crystallization is thought to be caused by a temperature gradient. After 20.0 ns, the crystallization slowly propagates, and crystallization almost stopped at 50.0 ns.

Fig. 1. Crystallization processes of amorphous Ge. (a) Atomic configurations and (b) trajectories ofatoms. The red lines show high-mobility atoms corresponding to a liquid state, and the green lines
show low-mobility atoms such as amorphous and crystal states. The temperature was set to be 1800
K and 1000 K at the top and bottom area, respectively.

2

JSST2022
Fig. 1(b) shows trajectories of atoms which are continuous lines obtained from the timedependent particle positions. The atoms in the liquid state have high mobility so that their trajectories fill the space. The trajectories of atoms with large mean square displacements areshown by red lines. At 14.0 ns, the large liquid region is locally formed. After 16.0 ns, thecrystallization proceeds toward the liquid layer and the volume of liquid decreases. This suggests that the liquid layer assists the rapid crystallization. The explosive crystallization observed experimentally is thought to occur in a similar process.
Fig. 2 shows the average temperature distribution in the MD cell. The
temperature of the right and left ends
was set to 1800 K and 1000 K, respectively, as indicated by a red broken line.
On the other hand, the actual temperature distribution is different. We heated
only the narrow area of 2 nm from far
right at 1800 K, but the high-temperature area with >1800 K expands up to
the half of the MD cell. This temperature distribution indicates that liquid
layers are easy to be formed near thehigh-temperature area. The liquid-me-

Fig. 2. The temperature distribution of the MD cell witha temperature gradient. For the simulation, the right andleft ends were set to 1800 K and 1000 K, respectively,
as indicated by a red line. The actual temperature is plotted as dots.

diated crystallization release heat andamorphous phase absorbed heat and melted. This react is self-sustaining propagation which
induce this broadened over 1800 K area.

4. Summary
The crystallization processes of amorphous Ge were examined by MD simulations using the
Tersoff interatomic potential in combination with the Langevin approach. Amorphous Ge washeat treated under a temperature gradient of 1800 K for the high temperature side and 1000 K
for the low temperature side. Rapid crystallization was observed under this heat treatment condition. A liquid layer was formed during crystallization, which was found to assist rapid crystallization. Explosive crystallization is expected to occur in a similar process.

3

JSST2022

References[1] G. C. Egan, Acta Mater 179 190 (2019).
[2] K. Kohno and M. Ishimaru, Jpn. J. Appl. Phys. 57, 095503 (2018).
[3] J. Tersoff, Phys. Rev. B 39, 5566 (1988).
[4] LAMMPS Molecular Dynamics Simulator. https://lammps.sandia.gov/.
[5] A. Stukowski, Modelling Simul. Mater. Sci. Eng. 18, 015012 (2010).
[6] OVITO User Manual – OVITO – Open Visualization Tool. https://www.ovito.org/d
ocs/current/.

4

JSST 2022

酸化タンタルを用いた ReRAM 素子の抵抗変化特性に対する電圧印加速度の影響
Effect of voltage application rate on resistance change
characteristics in Ta2O5-based ReRAM cells
石津 岳 1,2*, 西 佑介 1,
Gaku Ishizu1,2*, Yusuke Nishi11

舞鶴工業高等専門学校 電子制御工学科
豊橋技術科学大学 工学部 電気・電子情報工学科
1
Department of Control Engineering, National Institute of Technology, Maizuru College, Japan2
Department of Electrical and Electronic Information Engineering,
Toyohashi University of Technology, Japan*ishizu.gaku.ge@tut.jp
2

Abstract. In recent years, resistance change memory (ReRAM) has been expected as an alternative to SSD and HDD. In this study, resistance change characteristics in a Ta2O5-based ReRAM cell were investigated using a coupled simulation model in which a sinusoidal voltage isapplied to the cell with transport of diluted species, electric currents, and heat transfer solvers
in COMSOL Multiphysics. The combination of different voltage application rates was usedin the coupled simulation. The results indicated that the resistance change characteristics varied
according to positive and negative voltage application ratios. In other words, the ratio is an
important factor in controlling resistance change characteristics in the asymmetric Ta2O5-based
ReRAM cell.
Keywords: resistance change memory, oxygen vacancy, diffusion equation

1. まえがき
抵抗変化型メモリは、酸化物を金属で挟んだ単純な構造を有する．初めに抵抗変化素子に電圧印加を行うことで酸化物層内に酸素空孔という正の電荷を持つ格子欠陥により構成される導電性フィラメントが形成される．その後に電圧印加を繰り返す
ことで導電性フィラメントが断裂および再形成，すなわち素子の電気抵抗の変化を繰り返すため，不揮発性メモリとして動作する．酸化物としては主に二元系酸化物が多く，中でも Ta2O5（酸化タンタル）を用いた構造はパナソニック社が採用している．
Ta2O5 系の熱力学的安定性が理由とされているが，詳細は明らかにはされていない．

5

JSST 2022
本研究の目的は，Ta2O5 を抵抗変化層に，より酸素組成が小さい TaOx を酸素空孔貯蔵層に用いた素子のアナログ抵抗変化の駆動力を，電磁気・伝熱・希釈種輸送の連成シミュレーションを用いて調べることである．酸素空孔による導電性フィラメントの形状や導電率などを定量化し，抵抗変化特性の改善指針を探索することを目指す．

2. 研究背景・方法
本研究で扱うシミュレーションモデルでは，酸素空孔の移動に伴う 2 つの物理現象を前提とする．電圧の大きさに比例するドリフトと，酸素空孔の濃度勾配や温度に強く 依 存 す る 拡 散 で あ る ． 以 前 の 研 究 で は 、 白 金 (Pt) を 上 下 の 電 極 に 用 い た Pt/
TaOx/Ta2O5/Pt の積層構造を有する素子が示す独自の抵抗変化特性[1,2]を，このモデルを用いたシミュレーションで再現することに成功した[3,4]．負の電圧最大値を増大させると，抵抗変化特性の極性変化が起こる点に独自性がある．
本研究では，このモデルで Fig. 1 のように印加速度 0.5 V/s で電圧最大値を 1.5 V に固定し電圧印加を繰り返しても，抵抗変化特性の極性が変化しないことを確認した．
そこで，ドリフトと拡散との関係が抵抗変化層内の酸素空孔濃度(cv)の分布 (導電性フィラメントの形状) に与える影響を考察すべく，アナログ抵抗変化特性に対する電圧印加速度依存を調べた．

Fig. 1 Resistance change characteristics of multiple applications of the same voltage

3. 研究方法
酸素空孔濃度の拡散モデルは，以下の式で表される．
𝐸𝑎
𝜕𝑐𝑣
𝑍𝑞
= 𝐷0 𝑒 −𝑘𝑇 ∇ ∙ (∇𝑐𝑣 −
𝐸𝑐 )
𝜕𝑡
𝑘𝑇 𝑣
右辺の第 1, 2 項は，それぞれ拡散とドリフトを示す．ここで，𝐸𝑎 は活性化エネルギ
ー，𝑘はボルツマン定数，𝑇は温度，𝑍は酸素空孔イオン価数，𝑞は素電荷，𝐸は電界である．これら 2 項のバランスにより Ta2O5 中の𝑐𝑣 の分布が初期状態(導電性フィラメ
ントの形成状態)から変化する．この拡散方程式と，電流の連続式およびフーリエの
熱伝導方程式を，有限要素法をベースとした COMSOL Multiphysics を用いて連成さ
せ，電圧印加に伴う抵抗変化時の素子領域内の𝑐𝑣 , 𝑇, 𝐸各分布結果を得た．
これらのモデルが示唆するのは，印加電圧の時間変化が抵抗変化特性に大きな影響を及ぼすということである．そこで，本研究ではまず，電圧最大値を 1.5 V で固定したまま，正負の電圧印加速度を変化させた．

6

JSST 2022

4. 結果・考察
まず，正負共に電圧印加速度 R を同じだけ変化させた．この変化では抵抗変化特性の極性変化は見られなかった．しかし，電圧印加が終了した際の導電性フィラメント
近傍の濃度を示した Fig. 2 より，電圧印加速度が低下するほど濃度の拡散が進んでいないことが確認できた．

Fig. 2 Oxygen vacancy concentration when the same positive and negative voltage application
rates are different

また，負の電圧印加速度を 0.5 V/s に固定し，正の電圧印加速度のみを変化させた
際の電流電圧（I-V）特性を Fig. 3 に示す．正の印加速度＜負の印加速度では極性変化
が，正の印加速度＞負の印加速度では抵抗変化極性の変化はなく抵抗変化比の増大が起こった．一方，正の電圧印加速度を 0.5 V/s に固定し，負の電圧印加速度のみを変化させた結果である Fig. 4 からも同様に，正の印加速度＜負の印加速度では極性変化，正の印加速度＞負の印加速度では抵抗変化比の増大のみが起こることが確認できた．

Fig. 3 I-V characteristics under different positive voltage application rates

Fig. 4 I-V characteristics under different negative voltage application rates

7

JSST 2022
以上の結果を考察する．電圧印加速度の低下による導電性フィラメント内の酸素空孔濃度の増加は，電圧印加時間の延長によるドリフト効果の拡大と，発熱量増加で
TaOx 層からの酸素空孔供給量の増大が，導電性フィラメントからの酸素空孔の拡散量を上回ったものと考えられる．また，電圧印加によるドリフトと拡散とのバランスが同程度であれば，変化させる電圧印加速度の正負に関係なく抵抗変化特性の極性や抵抗変化比の傾向が等しくなることがわかった．つまり，酸素空孔供給層を含む本モ
デルでは，正と負のドリフト比を適切に調節することで抵抗変化特性を制御できるといえる．以上のことから，抵抗変化特性と導電性フィラメントの形状には明らかに強い相関があることが，シミュレーションの導入により確認することができた．

5. おわりに形成箇所が定まらない数 nm 程度の系を有する導電性フィラメントに連成シミュ
レーションを適用させることにより，実験では確認できない多くの知見が得られた．
抵抗変化特性の極性変化の程度とは無関係に，正負の電圧印加量の差を拡大させることが，抵抗変化特性を明確にする上で重要である．アナログ抵抗変化特性を改善する
ためには，導電性フィラメント内の酸素空孔濃度の増加による安定化の観点では，正負ともに電圧印加速度を低下させることが最も効果的であると考えられる．

参考文献[1] T. Miyatani, Y.Nishi, T. Kimoto: Two modes of bipolar resistive switching characteristics
in asymmetric TaOx-based ReRAM cells, MRS Adv. 4 (2019) 2601-2607.
[2] T. Miyatani, Y.Nishi, T. Kimoto: Methods of controlling operation modes in
Pt/TaOx/Ta2O5/Pt resistive switching cells, Jpn. J. Appl. Phys. 61 (2022) SM1006 1-6.
[3] Y. Nakamura, Y. Nishi: Theoretical Analyses of Analog Resistance Change Phenomena
in Pt/TaOx/Ta2O5/Pt cells, MEMRISYS 2021, P-11.
[4] Y. Nakamura, Y. Nishi: Simulation of analog resistance change characteristics in
Pt/TaOx/Ta2O5/Pt cells, Jpn. J. Appl. Phys. 61 (2022) SM1012 1-7.

8

JSST2022

Simulation of Panic Buying in a Disaster
Focusing on Consumer Information Transfer
Taku Ideguchi1*, Kazuo Furuta1,2, Taro Kanno1
Department of Systems Innovation, School of Engineering, The University of Tokyo
Resilience Engineering Research Center, School of Engineering, The University of Tokyo*ideguchi-taku@g.ecc.u-tokyo.ac.jp
12

Abstract. Panic buying in a disaster makes it difficult for many people to obtain supplies andcauses supply chain disruption. These days, due to the development of social networking services, the information transmitted by consumers often has a significant impact on panic buying.
Our previous research has analyzed the effect of consumer panic buying on the supply chain
by agent-based simulation, but it had some problems with the information transfer and it did
not take into account measures taken by consumers. Therefore, the objective of this study is todevelop a more realistic model and formulate measures that can be taken by consumers. The
results showed that transmitting counter-evidence from the hub agent in 15 or 20 days after theoccurrence of rumors can be the best measure to prevent supply chain disruption and maximize
consumer satisfaction. This simulation results could be used to help people or organizations who want to correct misinformation and reduce panic buying in a disaster.
Keywords: panic buying, consumer behavior, information transfer, agent-based simulation

1. Introduction
When a disaster occurs, in addition to the direct damage, secondary damage can be seen.
Hoarding of food and daily necessities caused by the uncertainty due to the disaster is a typicalexample of secondary damage. This hoarding behavior is called panic buying. These days, because of the development of social networking services, many consumers can disseminate information, which contains inaccurate information. The diffusion of incorrect information canlead to panic buying. In times of disaster when panic buying occurs, it is necessary to prevent
supply chain disruptions and ensure consumers have access to supplies.
From these backgrounds, Kawanishi extended Dulam’s model of panic buying [1] and did anagent-based simulation on supply chain disruption caused by panic buying of bottled water
during a disaster and proposed countermeasures that can be taken by the supply chain agents(SCAs) (hereinafter referred to as the previous study). In the previous study, he concluded thatthe impact of panic buying could be effectively mitigated by combining sales restrictions withemergency safety stock, sharing inventory information among SCAs, and provision of retailstore inventory information to consumers. However, this study has some problems. For

9

JSST2022
example, the previous study did not consider the measures that can be taken by consumers.
Therefore, in this study, a more realistic model was developed by extending the previous study,
to formulate realistic measures that consumers can take.

2. Method
2.1 Overview of the model
The product targeted in the model is bottled water. The model is based on the previous study
and consists of two parts: the supply chain model and the consumer model. The supply chain
includes three tiers: manufacturer, distributor, and retailer. The retailer sells the product to consumers. Additionally, there are two types of retailers: independent retailers and chain retailers.
Chain retailers contain the chain head and stores. The chain head buys products from the distributor and stores the product. Chain stores wholesale products to the chain stores. Each consumer agent represents a household. Consumer agents consume water every day and take purchase action when their stock runs low or runs out. In addition to this normal purchase action,
when a disaster occurs, consumers may panic buy. The conditions for panic buying are determined by four factors: the proportion of retail outlets restricting purchases, the amount of media
coverage, the number of rumors, and the proportion of consumers in the neighborhood engaging in panic buying behavior. To measure from both supply chain and consumers, measure thechange of the average stock in each tier of the supply chain and the change in consumers’satisfaction.
2.2 Model improvements
Some points in the consumer model used in the previous study are very simple. First, it isunrealistic to assume that rumors are spread uniformly to all consumer agents. Also, although
only one rumor is considered as triggering panic purchasing behavior, in reality, information
counteracting it is likely to spread. Another problem is that there is no mention of measuresthat can be taken by consumers. To address these problems, a scale-free network was set up as
a rumor diffusion network among consumers. This network was generated by using the Barabasi-Albert model [2]. In addition, as a minor change, the decision-making to engage in panic
buying behavior was changed to make probabilistically based on the four factors mentioned
above. Furthermore, the effect that the psychological impact of purchase restrictions decays
overtime was implemented.
2.3 Measures
Rumors in the previous study are assumed to be false information, and in addition to thisrumor, counter-evidence is implemented in the model. Counter-evidence is the information tocounteract false rumors, so this can be a countermeasure to rumors. It is assumed that rumors
have a facilitating effect and counterevidence has an inhibiting effect on panic buying. Counterevidence spreads through the information diffusion network among consumers in the same way
as rumors. Rumors were assumed to be false information about the availability of bottled water
during a disaster, and counter-evidence was assumed to be its correction.
Each consumer agent has one of three states of believing a rumor: easily believe the rumor,
hardly believe the rumor and neither. The opposite is true for the counter-evidence. Consumers
will know about a rumor or counter-evidence if there is more than a certain percentage of10

JSST2022
consumers around them who are transmitting rumor or counter-evidence, and they will decide
whether or not to believe it according to their characteristics.
Moreover, the reliability of counter-evidence is set on three levels: high, middle, and low. Thisparameter represents the effect of the source and the channel of the information and so on. The
reliability has an effect on the aforementioned criteria for deciding whether to know about arumor and counter-evidence. The rumor and counter-evidence assumed in this study are false
information about the availability of drinking water in the disaster and its correction. So, counterevidence was assumed to be more reliable than rumors, even if they have a low confidence
level.
Two measures were taken regarding the counter-evidence. The first measure was selecting
agents to initiate counter-evidence. The second measure was to find the optimal timing of initiate counter-evidence.
2.4 Simulation settings
To evaluate the effectiveness of the proposed methods, simulations were conducted for several
scenarios. The initiator of information was selected at random in the base scenario, while the
hub agent was the initiator in the other. The hub agent here means the agent who has the largest
number of links among surrounding agents. As for the timing of initiating counter-evidence, 5scenarios were considered: 1, 5, 10, 15, and 20 days after the rumor occurrence. The characteristics of consumer agents, and how likely they believe rumors or counter-evidence, are allocated in 2:6:2(easy to believe: neither: hard to believe). The simulation was run for 100
timesteps and one timestep represents one day. The disaster occurrence time was set to 70thtimestep. Measures were taken one step after the disaster occurrence. The two performance
metrics used are the inventory ratio and the consumers’ satisfaction. The former can be measured by reviewing the changes in inventory ratio in each supply tier. If the ratio is more thanone, the supply tier has sufficient inventory, and if the ratio is less than one, the supply tier isconsidered to be stock out. If multiple supply tiers are continuously stocked out, the supply
chain is regarded as disrupted. The latter can be measured by reviewing the changes in thenumber of satisfied consumers. The satisfaction of consumers is divided into the five categories: fully satisfied, highly satisfied, satisfied, lowlily satisfied, and not satisfied at all.

3. Results
Figure 1 shows the changes in the number of consumers in each satisfaction level in the selection of the agent who initiates counter-evidence. The result shows that initiation of the counter-evidence from the hub agent can be an effective measure to maximize the number of fully
satisfied consumers.
From the hub agent

Number of consumers

From random agent

full
high
middlelow
zerotimesteps

timesteps

Figure 1: Consumer satisfaction(random vs. hub)11

JSST2022
Figure 2,3 shows the changes in the consumers’ satisfaction and inventory sales ratio in thetiming of the initiation of counter-evidence. In figure3, red areas mean supply chain disruption.
This result shows that sending the counter-evidence one step after the rumor occurrence is thebest measure to minimize the effect of panic buying. However, in this simulation, the counterevidence is assumed to be more reliable than the rumor. If the rumor was an outright lie, it iseasy to make highly reliable counter-evidence, but in any case in the real world, rumors are
difficult to affirm to be false in such a short time. So, to deal with this problem, 4 measures
were considered excluding measures to disseminate counter-evidence immediately after the
rumor occurrence. As a result, sending counter-evidence 10 days after the rumor occurrence isthe best measure in terms of satisfaction, and 15 days are the best measure in terms of inventory
ratio.
10 steps

15 steps

Number of consumers

1 stepfull
high
middlelow
zerotimesteps

timesteps

timesteps

Figure 3: Consumer satisfaction (change the timing of counter-evidence initiation)
10 steps

15 steps

Inventory sales ratio

1 stepretailerchain head
distributor
manufacturer

timesteps

timesteps

timesteps

Figure 2: Inventory sales ratio of SCAs

4. Conclusion
This research mainly focused on the effect of counter-evidence. When the reliability of counter-evidence is high, it was confirmed that an effective countermeasure is to start diffusing
counter-evidence from consumer agents who have a strong influence on their surroundings. It
can also be said that if counter-evidence is more reliable than rumors can be disseminated,issuing counter-evidence immediately after the occurrence of a rumor can be a realistic andeffective measure against panic buying. If counter-evidence is not reliable as a rumor, it is goodto spread counter-evidence in 10 or 15 days after the rumor starts spreading. This simulation
results could be used to help people or organizations who want to correct misinformation and reduce panic buying in a disaster. In the future, implementing the dynamics of consumer attitudes toward information could improve the model.

References[1] R. Dulam, K. Furuta, T. Kanno: Development of an agent-based model for the
analysis of the effect of consumer panic buying on supply chain disruption due toa disaster, Journal of Advanced Simulation in Science and Engineering, vol.7(2020), Issue 1, 102-116.
[2] R. Albert, H. Jeong and A. L. Barabasi: Diameter of the world-wide web, Nature,
vol.401(1999), 130-131

12JSST 2022

脳神経細胞群のマルチスケール
同期シミュレーション法の開発
Development of Multiscale Synchronic Simulation
Method for Brain Neuronal Groups
池田 翔 1*
Sho Ikeda1*

1

同志社大学大学院 生命医科学研究科 医工学・医情報学専攻
1

Department of Biomedical Engineering and Informatics,

Graduate School of Life and Medical Sciences, Doshisha University*ctug0008@mail4.doshisha.ac.jp

Abstract. Synchronization of neuronal firing in the brain has attracted much attention in recent years. In this study, we aimed to develop a mathematical model to reproduce synchronizations using a neuronal membrane potential model as microscale and a micro model andcellular automata structure as a macroscale for multiscale simulation. The micro model wasdevised by averaging the membrane potentials of all cells on the microscale and modeling theobtained pulse waveform in Fourier series. In the results of microscale and macroscale simulation, I found synchronic frequency of each scale consistent. Therefore, multiscale synchronic simulations were obtained from this study.
Keywords: Brain simulation, multiscale, synchronization model, cellular automata

1. 緒言
脳内で生じる神経細胞の発火の同期現象は，脳活動の分析等に役立ち，近年注目されている[1]．先行研究で脳の同期シミュレーションがある[2]が，計算コスト削減に向けたモデルの簡易化や数値計算法として正しいと言えるものが少ないことが分かった．そこで本研究では，Fig.1 の概要図のように発火の同期現象を再現する数理モデ
ルとして，ミクロスケールでは神経細胞間の生理学的現象であるシナプス電流とSTDP のモデルを取り入れた神経細胞の膜電位モデルを構築し，マクロスケールでは，
ミクロスケールで得られた同期現象をモデル化したミクロモデルとセルオートマト
ンを用いることによる，マルチスケール同期シミュレーション法の開発を目的とした．

13

JSST 2022

Fig.1

Outline figure of this study

2. 数理モデル
ミクロスケールで扱うモデルを式(1)~(4)に示す．式(1)(2)を膜電位モデル，式(3)をシナプス電流モデル，式(4)をシナプス荷重モデルとした．
𝑖𝑖𝑣̇ = 0.04𝑣𝑖 2 + 0.5𝑣𝑖 + 140 − 𝑢𝑖 + 𝐼𝐷𝐶− 𝐼𝑠𝑦𝑛+ 𝑠𝑤𝑖{ 𝑖𝑢̇𝑖 = 𝑎𝑖 (𝑏𝑖 𝑣𝑖 − 𝑢𝑖 )𝑣 ← 𝑐𝑖𝑢𝑖 ← 𝑢𝑖 + 𝑑𝑖𝑖𝑓 𝑣𝑖 ≥ 30[𝑚𝑉] 𝑡ℎ𝑒𝑛 { 𝑖

(1)
(2)

𝑁

𝑖
𝐼𝑠𝑦𝑛
=

𝐽
𝑗
∑ 𝑠𝑗 (𝑡)(𝑣𝑖 − 𝑉𝑠𝑦𝑛 )
𝑁−1

(3)

𝑗(≠𝑖)
𝑁

𝜂
𝑠𝑤𝑖 =
∑ 𝑆𝑇𝐷𝑃𝑗𝑋,𝑌 (𝑡𝑝𝑜𝑠𝑡 − 𝑡𝑝𝑟𝑒 )
𝑁−1

(4)

𝑗(≠𝑖)

式(1)において，𝑣̇ 𝑖 ，𝑢̇ 𝑖 を時間𝑡[msec]に対する𝑖番目のシナプス後細胞の膜電位𝑣𝑖 [mV]の変化量，リカバリー変数𝑢𝑖 の変化量とした[2]．シナプス後細胞𝑖への外部入力項を，𝑖
𝑖
視床入力𝐼𝐷𝐶，シナプス電流𝐼𝑠𝑦𝑛，シナプス荷重𝑠𝑤𝑖 とした．式(2)で𝑣𝑖 ≥ 30[mV]のとき，シナプス後細胞𝑖の発火を意味する．パラメータ𝑎𝑖 , 𝑏𝑖 , 𝑐𝑖 , 𝑑𝑖 で発火形状が指定できる[2]．
𝑖
シナプス電流𝐼𝑠𝑦𝑛
は，シナプス前細胞𝑗が発火することで，シナプス後細胞𝑖へ流れる電流である．これは軸索の長さに依存して伝達遅延が生じる．𝑁を全細胞数，𝐽をシナプス間の結合強度，𝑠𝑗 (𝑡)をシナプス前細胞𝑗の発火により放出される神経伝達物質
𝑗がシナプス後細胞𝑖の樹状突起受容部へ結合する確率，𝑉𝑠𝑦𝑛 をシナプス前細胞𝑗の逆転電位とした[3]．シナプス荷重は，神経細胞間で生じる STDP を膜電位の変化に反映させる役割を担う．STDP は，シナプス前後細胞の細胞種類と発火時間に依存し，細胞間のシナプスの強度を変化，保持する性質を持ち，学習や記憶，認知の機能の役割を担う．𝑡𝑝𝑟𝑒 , 𝑡𝑝𝑜𝑠𝑡 をシナプス前，後細胞の発火時刻，𝑋, 𝑌をシナプス前，後細胞の種類，𝜂を学習倍率，𝑆𝑇𝐷𝑃𝑗𝑋,𝑌 (𝑡𝑝𝑜𝑠𝑡 − 𝑡𝑝𝑟𝑒 )を𝑋𝑌間の STDP 関数とした[4]．
ミクロスケールシミュレーション結果で得られた周期的な同期現象時において，全細胞の膜電位を平均化し，得られた同期波形に対して，フーリエ級数を用い，ミクロモデルを構築した．マクロスケールにおけるセルオートマトン構造では，１つのセル

14JSST 2022を１つのミクロモデルとし，隣接するセル間では同期発生によるミクロモデルの位相シフト効果を取り入れた．以下の式(5)にミクロモデル，式(6)にシフトモデルを示す．
3

𝑉𝑖 (𝑡) = 𝑎0 + ∑[𝑎𝑛 cos{𝑛(𝜔𝑡 − 𝜑𝑖 )} + 𝑏𝑛 sin{𝑛(𝜔𝑡 − 𝜑𝑖 )}]

(5)

𝑛=1
𝑁𝑖

𝐾
𝜑𝑖 =
∑ 𝐻𝑚 (𝑉𝑗𝑚 − 𝑉𝑖 )
𝑁𝑖

(6)

𝑚=1

式(5)において𝑉𝑖 (𝑡)を𝑖番目のセルにおける時間𝑡[msec]の膜電位[mV]，𝑎0 , 𝑎𝑛 , 𝑏𝑛 , 𝜔を関数パラメータとした．式(6)の𝐾を結合強度，𝑁𝑖 を𝑖番目のセルと繋がるセル数，𝑚を隣接するセルの座標番号，𝑉𝑗𝑚 を𝑚番目の隣接するセル𝑗の膜電位とした．同期判定𝐻𝑚を同期発生有(𝐻𝑚 = 1)，同期発生無(𝐻𝑚 = 0)とした．

3. シミュレーション結果と考察

Fig.2

A result of microscale simulation

Fig.4

A result of time-frequency

Fig.6

Number of synchronous firing cells

Fig.3 Macro model

Fig.5 A result of macroscale simulation

Fig.7 A result of time-frequency (𝐾 = 1.0)

ミクロスケールシミュレーションの条件として，興奮性細胞と抑制性細胞の比率

15JSST 2022𝑖
4:1 とした全細胞数𝑁 = 100，視床入力𝐼𝐷𝐶= 4.0，シナプス間の結合強度𝐽 = 10，学習倍率𝜂 = 2.0とした．シミュレーション結果を Fig.2 に示す．Fig.2A の点プロットは，発火した神経細胞を意味する．興奮性細胞と抑制性細胞で同期現象の発生有無を同期判定手法で確認し，その結果を Fig.2B に示す．Fig.2B の点プロットは，同期発生を意味する．全細胞の膜電位を平均化し，フーリエ級数を用いてミクロモデルを構築した．構築した式(5)のミクロモデルにおいて，Fig.3 に時間に対する膜電位𝑉𝑖 (𝑡)を示す．
Fig.3 において，１周期中の最大ピーク値が Fig.2B の点プロットに相当する．さらにウェーブレット解析を用いて時間周波数解析を行った．その結果を Fig.4 に示す．Fig.4
から一定周期かつ 2 つの周波数を持つ波があることが分かる．ここから，同期発生を示す波において，同期周波数が約 10[Hz]を示していることが分かる．
マクロスケールシミュレーションの条件として，100×100×100 の 3 次元セルオート
マトン構造，初期時間𝑡[msec]を[0.1,100]の一様乱数，結合強度𝐾 = 0.5, 1.0, 1.5とした．
Fig.5 にシミュレーション結果，Fig.6 に時間に対する同期発火セル数を示す．Fig.6 の結果に対して，ウェーブレット解析を用いて時間周波数解析を行った．その結果をFig.7 に示す．Fig.7 から初期では複数の周波数を持つ特徴が見られたが，十分時間が経過すると，約 10[Hz]に収束した．Fig.4 と Fig.6 の同期周波数を比較すると，ミクロ，
マクロスケール同期シミュレーションの同期周波数が一致するため，本研究で開発したシミュレーション法により，マルチスケール同期シミュレーションが再現できた．
ここで，この Fig.7 の収束の原因に着目する．式(6)のシフトモデルによる位相シフ
トは，𝑖番目のセルに対して，隣接するセルの同期発火発生と𝑖番目のセルとの膜電位の差分で決まり，膜電位の差分が大きい場合，位相シフト量も大きくなる．ここで初期値が一様乱数であることから式(5)の膜電位𝑉𝑖 (𝑡)も一様乱数であり位相シフトが生じる．これにより，膜電位𝑉𝑖 (𝑡)の位相は，同期発火した隣接セルの膜電位に近い値へシフトするため，同期発火セル数が 1 つの周期に収束していくと考える．

4. 結言
本研究で開発したシミュレーション法により，マルチスケール同期シミュレーションが再現できた．今後は，結合強度𝐾と同期発火セル数の変化の関係性を調査し，これを基に 3 次元空間内に 2 種類の結合強度𝐾を取り入れ，同期発火パターンの変化を調べていく．

参考文献[1] S. N. Kbah: A Computational Model of the Brain Cortex and Its Synchronization, BioMed Research International, (2020).
[2] E.M. Izhikevich: Polychronization: Computation with Spikes, Neural Computation,
vol:18, no.2 (2006), 245-282.
[3] S.Y. Kim, Coupling-Induced Population Synchronization in An Excitatory Population ofSubthreshold Izhikevich Neurons, Cognitive Neurodynamics, vol:7 (2013), 495–503.
[4] S.Y. Kim, W. Lim, Effect of Interpopulation Spike-Timing-Dependent Plasticity on
Synchronized Rhythms in Neuronal Networks with Inhibitory and Excitatory Populations, Cognitive Neurodynamics, vol:14, no.4 (2020),p. 535-567.

16

JSST 2022

Numerical prediction of water-ice
phase-change with natural convection
including density inversion region
Kazuki Nishimoto1,* , Daisuke Toriu2 , Satoru Ushijima2
1 Graduate School of Engineering, Kyoto University
2 ACCMS, Kyoto University∗ nishimoto.kazuki.35n@st.kyoto-u.ac.jp

The full paper of this paper has been published in Transaction of the Japan Society for
Simulation Technology.
DOI: https://doi.org/10.11308/tjsst.15.27

17

JSST 2022

フラクタル的に成長する埃の堆積過程の考察
Consideration of Dust Accumulation Process ofGrows Fractal-wise
大谷 昌生 1*, 金井 徳兼 1
Masaki Ohtani1*, Norikane kanai 1

1
1

神奈川工科大学 大学院 工学研究科 電気電子工学専攻

Department of Electrical and Electronic Engineering, Graduate School of Engineering,
Kanagawa Institute of Technology*ootani2014@

ele.kanagawa-it.ac.jpu

Abstract. Natural landscapes and existing shapes can be treated with the idea of fractals, which
do not have a characteristic length and can be considered by self-similarity. Dust is also considered to have a self-similar fractal shape without a characteristic shape. Thought that modeling and simulating randomly growing dust will lead to non-contact measurement of the threedimensional accumulation amount of dust. In this paper, image generation generating and analyzing of images of 3D dust deposit is reported using simulations.The correlation between
fractal dimension and spatial occupancy is used to examine the effectiveness of the 3D dust
deposition model.
Keywords: Dust, Fractal, Image analysis

1. はじめに
自然界の風景や存在する形は、特徴的な長さを持たず自己相似性で考えることができるフラクタルの考え方で取り扱うことができる。特徴的な長さを持たない形の例は、積乱雲やコッホ曲線のほかにもたくさん知られている。例えば、海岸線や山の起伏や川の形などが挙げられる[1]。人工的な構造物、鉄塔の骨組みがフラクタル形状を示す
現実もある。埃もまた生活空間内に発生すると小さな埃同士が集合し大きな埃へと成長する。この成長過程は、水蒸気が集まってできる雲の立体形状に酷似していると考え、大きな埃の形成もフラクタルであると著者らは考える。著者らは、埃の堆積量や３次元的な空間占有率を非接触的に計測する方法について取り組んでいる。ランダムに成長する埃をモデル化しシミュレーションすることで、埃の３次元的な堆積量の非接触的な計測に繋がると考える。

21

JSST 2022
本論文では、埃の特性を観察し、ランダムに成長する埃の堆積過程をモデル化する。
微小な埃の集合過程を 2 次元平面でランダムに成長させ積層化した 3 次元的に成長する埃堆積モデルを提案する。ランダムに成長する 3 次元埃堆積過程をシミュレーシ
ョンし、３次元的な埃堆積画像を示す。シミュレーションから得た 3 次元埃堆積画像を解析し、フラクタル次元と空間占有率の相関から、3 次元埃堆積モデルの有効性を考察すると共に、埃がフラクタルな成長過程であるかを議論する。

2. 埃の性質とフラクタル
埃の構成成分から、衣服などから抜け落ちるポリエステルなどの合成繊維が主である。ポリエステルなどの合成繊維は、微量に帯電しており、埃の同士に引力が発生しながら空気中を浮流する。一般住宅では、埃の帯電により埃同士が静電気力により引き合いながら床に堆積すると考える。また、繊維質形状で微小な埃同士が接触すると絡み合う。時間とともに埃の堆積量が増加する過程は、3 次元にランダムに成長する図形であると著者らは考える。
フラクタルとは、海岸線や雲などの特徴的な長さを持たない図形や構造、現象などの総称である。画像処理におけるフラクタル次元解析では、画像濃度曲面の起伏の複雑さを表す特徴量となる。画像処理によるフラクタル次元解析を活用することで、構造物や製品の表面を数値化し、品質評価に活用されている[2]。解析方法は、デジタル
画像の RGB 値やグレースケール下での濃度レベルをピクセル毎に数値化、解析範囲を設定し、ボックスカウンティング法やピクセルカウンティング法[3]を用いてフラクタル次元の解析できる。埃の堆積画像から画像処理によるフラクタル次元により埃
堆積情報を検出できると考える。

3. 埃の堆積シミュレーションの作成
2 次元平面にランダムに生成する埃モデルを考える。x-y の 2 次元平面にランダムに現れる埃の成長開始点を点 𝑃" (𝑥0 , 𝑦0 )とする。次点で発生する埃もまたランダムに発生することから、次点座標 𝑃) (𝑥) , 𝑦) )を示す式はランダム関数により、
𝑥) = 𝑥" + 𝐾𝑎." ， 𝑦) = 𝑦" + 𝐾𝑏1"とする。ここで用いる𝑎.2 , 𝑏12 は、−0.5〜0.5[㎜]範囲のとし、ランダム変数𝐾は成長倍率定数 (0 < 𝐾)である。ここで示すモデルは、成長倍率定数𝐾を 1 とするが、シミュレーション時の観測範囲を変化させる時、成長倍率定数𝐾を調整することで観測範囲の拡大縮小ができる。更に、次々点にランダムに発生する埃の座標 𝑃4 (𝑥4 , 𝑦4 )は、
𝑥 = 𝑥) + 𝐾𝑎.) ， 𝑦4 = 𝑦) + 𝐾𝑏1)となる。また、𝑛 回ランダムに成長する埃の座標 𝑃6 (𝑥6 , 𝑦6 )は、
66𝑥6 = 𝐾 7 𝑎.2 ， 𝑦6 = 𝐾 7 𝑎12
28"

28"

22

JSST 2022となる。ランダムに現れた埃の座標を描画（プロット）すると図 1(a)のようになる。
ランダムに発生した座標を粗視化することで、埃の成長をデジタル化する。ここで粗視化する格子サイズを一辺 0.1[㎜]とし、後に画像化するピクセルサイズと同じにする。次にランダム発生する埃の座標の最大距離が 0.5[㎜]のため、隣接した埃同士は絡み合う。一方、離れた座標で発生した埃は孤立するが、後にその座標に隣接した埃が発生することで、成長することができる。
3 次元的な堆積モデルは、2 次元平面に 𝑛 回成長した埃モデルを 𝑁 層作成し、重ね合わせることで 3 次元的な埃の堆積成長モデルを作成する。図 1(b)に、3 次元的な埃
の堆積成長モデルを示す。合計の成長回数𝑁′は、
𝑁′ = 𝑛 × 𝑁となる。x-y 方向 𝑛 回成長と z 方向 𝑁 層の積で成長回数 𝑁 @ を決定する。
Ｎ
・
・
・
・
・
Ｎ=３
・
・
・
・
Ｎ＝２
・
・
・
・
・
Ｎ=１

y

・
・

!' #' , %'

・

・
z
x

!& #& , %&・・!" #" , %"!) #) , %*

・・
・
・ ・

y

x

!( #( , %(

(a)

(b)
図 1

埃堆積モデルのイメージ.

(a) 2 次元方向に成長する埃の成長過程, (b) 3 次元埃堆積モデルの積層イメージ.

4. 埃の堆積過程のシミュレーション
3 項で検討した 3 次元的にランダムに成長する埃堆積モデルのシミュレーション
し、埃の堆積画像を生成する。前述の通り、シミュレーションで成長した埃 1 ピクセ
ルサイズは 0.1[mm]に相当し、 𝑧 軸方向の埃の重なりを、グレースケールの濃度値で画像を生成する。1 層の成長回数 𝑛 を 4 万回毎、4 万〜20 万回変化させ 4 層作成し、積層化した埃堆積モデルを画像化した。5 種のシミュレーション結果の埃堆積画像を
図 2 に示す。
埃のランダム的な堆積過程シミュレーションの結果をから生成した、堆積画像(図2)の空間占有率−フラクタル次元と空間占有率−成長回数の相関を図 3 に示す。フラク
タル次元解析には、ボックスカウンティング法を用いた。ボックスカウンティング法
は、対象とする画像を 1 辺が𝑟の正方形で覆ったとき，その中に対象とする点や線が含まれる正方形の数 𝑁(C) を数え，𝑟を変化した時のlog(𝑟) に対する logG𝑁(C) H の傾きから，フラクタル次元𝐷を算出する方法である．図 3 より、提案する 3 次元埃堆積モデ
ルの埃の成長回数と空間占有率及びフラクタル次元の比例関係が確認できたことから、埃の堆積画像のフラクタル次元から堆積量の推定できると考える。

23

JSST 2022

図 2

シミュレーションから生成した埃の堆積画像サンプル.
（埃の成長回数 左から𝑛 =4, 8, 12, 16, 20[万回]）.

2.80

28

20

2.70

16
2.65
12
2.60

8

y = 0.4x - 4.80

2.55

成長回数[万回]フラクタル次元24

y = 0.004x + 2.53

2.75

4

2.50

0
0

10

20

30

40

50

60

70

空間占有率[％]
空間占有率−フラクタル次元

図35.

空間占有率−成長回数

3 次元埃堆積モデルの空間占有率とフラクタル次元及び成長回数の相関.

まとめ

本報告では、埃の特性から堆積過程がランダムに成長する 3 次元埃堆積モデルを提案した。ランダム関数を活用し、2 次元平面に成長する埃を積層化することで、３次元的な埃堆積モデルを検討した。3 次元的な埃堆積モデルを活用し、粗視化したデ
ジタルデータに変換、グレースケール画像をシミュレーションから生成し、空間占有率とフラクタル次元の解析した。埃の堆積画像の解析結果から、埃の堆積形状には、フラクタル次元と空間占有率の相関があり、提案した 3 次元埃堆積モデルが有効であることを示した。

参考文献[1] 高安秀樹：フラクタル新装版,朝倉書店,(2010).
[2] 中山寛, 曽根光男,高木幹雄：フラクタル次元と低次統計量とを用いた気象衛星
N0AA 画像の解析とその評価, 情報処理学会論文誌 30 (1), pp.91-100, (1989).
[3] 宇田川義夫：トンネル切羽画像のフラクタル解析による地質評価手法の開発,
応用地質学会誌, 第 48 巻, 第 3 号, pp.116-125 頁, (2007).

24

JSST 2022

第五世代移動通信システムにおける基地局の配置問題
Base station placement problem in
5th generation mobile communication systems
紅林 照也 田村 裕
Shoya Kurebayashi Hiroshi Tamura

中央大学大学院 理工学研究科 電気電子情報通信工学専攻

Chuo University Graduate School of Science and Engineering
Department of Electrical,Electronic,Information and Communication Engineering
Tamura laboratory
a18.xy4m@g.chuo-u.ac.jp

Abstract. Currently,research is being conducted toward the spread of 5th generation mobile
communication systems. However, since 5G uses a high frequency band, it is necessary toconsider the communication range and obstracles. Applying the museum problem, the optimal arrangement of the antenna is required by triangulation and coloring. So, the optimum
arrangement of the antenna is simulated by dividing the area and coloring. Therefore, obstacles are placed in the area and triangulation is performed at each vertex. The coloring was
performed in descending order of the degree of each vertex,the method of reduce the numberof specific coloring points, the method of area,and the method of randomly coloring. The
method to reduce the number of specific coloring points is able to arrange a suitable base station.
Keywords: Computer simulation, Graph coloring, 5G, arrangement of antennas

25

JSST 2022

1. はじめに
近年、第 5 世代移動通信システムの普及に向け研究が行われており、5G 対応エリアも拡大してきている。5G では高周波数帯を使用するため 4G と比べて通信範囲が小さくなり、直進性が高く障害物に遮断されやすくなってしまう。そのため、5G ア
ンテナの配置には通信範囲や障害物による通信の影響を考慮しながら、なるべく少ない数のアンテナでエリア全体を被覆できるように適当な配置場所を選ぶ必要がある。
[1]
本研究では、計算幾何学及びグラフ彩色理論に基づいて、平面における障害物エリアを作成し、三角形分割を行い、最適で最も少ないアンテナの配置数を検討する。

2. シミュレーション概要
高密度に効率よく基地局を配置するには、美術館問題[2]という問題が応用できる。
美術館問題とは以下の条件のもと監視カメラの台数をできるだけ少なくする問題である。
・平面上に穴のない n 角形
・この領域全体を見渡せるように監視カメラを設置する・監視カメラは 360°監視できる
平面上に穴のない n 角形は三角形分割を行うことで外平面グラフとなり、外平面グ
ラフの点は 3 色で塗分け可能である。この 3 色の中で一番少ない色の点で監視カメラを配置すればよい。しかし、美術館問題は領域全体を見渡せるのであれば、距離は問
題ないが、5G アンテナの配置では周波数の関係で通信距離の制限がある。そのため通信範囲の制限の下で美術館問題を応用する。
図 2.1 のいくつかの障害物を含む領域を多数の正六角形エリアを作成し、この六角形エリアにおいて各頂点で三角形分割を行う（図 2.2）
。次に、各頂点で赤、緑、青、黄、桃、黒色の順番で貪欲彩色を行い、赤、緑、青色のそれぞれの点で 5G アンテナを配置し、通信範囲を周囲の辺の長さを半径とする円で電波は障害物に遮断されるという条件の下での通信範囲を図示する(図 2.3)。この時、障害物エリア内で被覆されていない場所が存在するため、その場所に新たなアンテナを配置する充足化という操作を行うことで障害物エリア全体が被覆され、赤、緑、青色の中で最も少ない色数が適当なアンテナの配置数となる。
彩色の方法として、ランダム、各頂点において次数の大きい点から順に彩色する総頂点次数の方法、各頂点において次数の大きい順から青色だけで隣り合わないように彩色し、その後次数の小さい順から貪欲彩色をする青最小の方法、被覆する面積の大
きい順から彩色する方法の 4 つで行う。また、障害物エリアの周囲を赤、緑、青色で固定する場合としない場合、障害物エリアを図 2.4 のように内部を正六角形で分割しないで周囲を赤、緑、青色で固定する場合としない場合の 4 パターンでシミュレーシ
ョンを行い、アンテナの配置数が最適となるか検討する。

26

JSST 2022

図 2.1 障害物エリア

図 2.2 三角形分割

図 2.3 充足化

図 2.4 エリアの分割なし

3. シミュレーション結果
障害物エリアの周囲を赤、緑、青色で固定しない場合を①、する場合を②、複数の正六角形エリアに分割しないで周囲を赤、緑、青色で固定しない場合を③、する場合
を④としてそれぞれの最小色数の点数のシミュレーション結果を表 3.1 に示す。また、ランダムは 10 回シミュレーションを行い、最小と最大の点数の結果を示す。
表 3.1 シミュレーション結果

ランダム
総頂点
面積
青最小

①
87~95
96

②
81~84
84
87
81

③
76~81
81
84
73

④
76~82
79
80
73

①と②を比較すると最小色数の点数をかなり少なくすることができた。また、青最小の方法が最も少ない点数の彩色方法となった。③と④を比較すると最小色数の点数の差はほとんどない。また、こちらも青最小の方法が最も少ない点数の彩色方法となった。

27

JSST 2022

4. まとめ
5G アンテナの最適配置を求めるために、ランダム、総頂点次数の方法、青最小の方法、面積の大きい順の 4 つの彩色方法と障害物エリアを複数の正六角形エリアに分けて彩色を行う方法、分けないで彩色を行う方法、これらを使って最小色数を検討した。
青最小の方法は充足化前の青色の彩色点数をできるだけ少なくなるようにアルゴ
リズムを作成した。まず、次数の大きい順に青色だけで隣り合わないように彩色を行い、次に次数の小さい順に貪欲彩色を行う。この貪欲彩色の際、青色の点において青色より優先度の高い赤、緑が彩色できる場合には赤または緑色を彩色する。この操作を行うことで青色の点を減らしていくことができる。④の時の充足化前の青色の彩色
点数を表 4.1 に示す。
表 4.1 充足化前の青色の彩色点数
ランダム
総頂点
面積
青最小
青色の彩色数
69~71
81
74
53
これより充足化前の青色の彩色点数を少なくできていることがわかり、あらかじめ色数を少なくしておくことが有効だと考える。
次に③と④は最小色数の点数はともに 73 であり、充足化前の青色の彩色点数は③は 45、④は 53 であった。③の場合はさらに充足化前の点数を少なくできたが、充足化後の最小色数の点数が同じなためで③では充足化で新しく追加する点が多くなる。
また、①と②から障害物エリアの周囲を固定した場合のほうが最小色数の点数は少な
いため、障害物エリアの周囲を固定するのは有効だと考える。
以上の結果から、充足化前の彩色点数をあらかじめ少なくなるように彩色することと、周囲の色は固定して彩色を行うことで適当なアンテナの配置数が求められた。また、さらにアンテナの配置数を減らすために点同士の距離で工夫してみるなどで試し
てみたい。

参考文献
[1] 王暁飛:”第五世代移動通信システムにおけるアンテナ配置”、日本シミュレーション
学会多次元移動通信網研究会講演論文集、JSST-MDMC2020-1-3、2020 年 12 月
[2] ガスコン研究所：「美術館定理」
http://gascon.cocolog-nifty.com/blog/2006/06/post_1fcf.html

最終閲覧日 2022 年 6 月 9 日

28

JSST 2022

Algorithm for Detour Hinge Vertex Problem
of Circular-arc Graphs ∗Tomonari Izumi1 , Syoma Nameki2 , Yoko Nakajima1 , Hirotoshi Honma1,* , Toshihiro Fujito3
1 Department of Creative Engineering, National Institute of Technology Kushiro College
2 Department of Electronics and Information Engineering, Hokkaido University
3 Department of Computer Science and Engineering, Toyohashi University of Technology∗ honma@kushiro-ct.ac.jp

Abstract. Consider a graph G = (V, E) with vertex set V and edge set E. The distanceδG (x, y) is defined as the length of the shortest path between vertices x and y in G. A vertexu ∈ V is a hinge vertex if there exist two vertices x, y ∈ V − {u} such that δG−u (x, y) >δG (x, y). The detour degree of vertex u is defined det(u) = max{δG−u (x, y) | δG−u (x, y) >δG (x, y), x, y ∈ N(u)}. For a hinge vertex u, the detour degree of u is defined det(u) =max{δG−u (x, y) | δG−u (x, y) > δG (x, y), x, y ∈ N(u)}. The detour hinge vertex problem aims
to determine the hinge vertex u that maximizes det(u) in G. In this paper, we propose aneﬃcient algorithm for the detour hinge vertex problem on circular-arc graphs.
Keywords: design and analysis of algorithms, detour hinge vertex problem, intersection
graphs, circular-arc graphs

1. Introduction
Consider a simple undirected graph G = (V, E) with vertex set V and edge set E. (Throughout the paper, n is the number of vertices in a graph.) Let G − u be a subgraph induced by
the vertex set V − {u}. The distance δG (x, y) is defined as the length (i.e., the number ofedges) of the shortest path between vertices x and y in G. Chang et al. [1] defined u ∈ V tobe a hinge vertex if there exist two vertices x, y ∈ V − {u} such that δG−u (x, y) > δG (x, y).
Finding all hinge vertices of a given graph is called the hinge vertex problem. This problem
has applications to the improvement of stability and robustness in communication network
systems [2].
Let U be a set consisting of all hinge vertices of G. The neighborhood of u, denoted by N(u), is the set of all vertices adjacent to u. Honma et al. [3] defined det(u) =max{δG−u (x, y) | δG−u (x, y) > δG (x, y), x, y ∈ N(u)} for u ∈ U as the detour degree of u. Thatis, det(u) means the degree to which a path between x and y becomes longer upon removal∗This work was partially supported by JSPS KAKENHI Grant Number 19K11834, 20K11093, and Cooperative Education/Research Project between Toyohashi University of Technology and National Institute ofTechnology.

29

of a hinge vertex u from G. The detour hinge vertex problem aims to determine the hingevertex u that maximizes det(u) in G. This problem is motivated by practical applications,such as network stabilization under a limited cost [2].
In this paper, we propose an eﬃcient algorithm for the detour hinge vertex problem on
circular-arc graphs. Circular-arc graphs are useful in modeling periodic resource allocation
problems in operations research. They also have applications in diﬀerent fields such as
genetic research, traﬃc control, computer compiler design and statistics [4]. Thus, circulararc graphs have been studied extensively from both a theoretical and algorithmic point ofview. Every articulation vertex obviously has a maximum detour degree. Therefore, in thispaper, we assume that a circular-arc graph does not include any articulation vertices, i.e., isbiconnected.

2.

Definition

We first illustrate the circular-arc model before defining the circular-arc graph. Consider aunit circle C and a family A of n arcs A1 , A2 , . . . , An along the circumference of C. Eacharc Ai has two endpoints, left endpoint ai and right endpoint bi , respectively, and is denotedby Ai = [ai , bi ]. The left endpoint ai (resp., right endpoint bi ) is the last point of Ai that
we encounter when walking along Ai counterclockwise (resp., clockwise). Without lossof generality, the coordinates of all left and right endpoints are distinct and are assigned
clockwise with consecutive integer values 1, 2, . . . , 2n. The arc numbers i, j are assigned
to each circular-arc in increasing order of their right endpoints bi ’s, i.e., Ai < A j if bi <b j . Note that arc Ai with ai > bi is called a feedback arc. The geometric representation
described above is called a circular-arc model. Figure 1(a) illustrates a circular-arc modelCM, consisting of 12 arcs (note that A1 and A2 are feedback arcs). Table 1 shows details ofCM.
A graph G = (V, E) is called a circular-arc graph if there exists a family of circular-arcs
A = {A1 , A2 , . . . , An } such that there is a one-to-one correspondence between vertex i ∈ V
and the circular-arc Ai ∈ A in such a way that an edge (i, j) ∈ E if and only if Ai intersects
with A j in CM. Figure 1(b) illustrates the circular-arc graph G corresponding to CM shownin (a). In this example, all hinge vertices of G are vertices 1, 2, 5 and 11. The detour degreeof each hinge vertex is det(1) = det(2) = 4, det(11) = 3, and det(5) = 5, respectively. If
hinge vertex 5 is removed from G, the distance between vertices 3 and 7 lengthens from 2
to 5. Thus, the detour degree of vertex 5 is det(5) = 5, and vertex 5 is a maximum detour
hinge vertex in G. We next introduce an extended circular-arc model(ECM) constructedfrom a CM. We first cut CM at point between 1 and 2n, next unroll onto the real horizontal
line. Each arc Ai = [ai , bi ] in CM is also changed to horizontal line segment. For each Ai ,1 ⩽ i ⩽ n, copies Ai+n are created by shifting 2n to the right. Honma et al. [6] presented an
O(n) time algorithm for the hinge vertex problem on circular-arc graphs by using ECM.

30

A1
A11

A12
24

2

12

4

2

4

C

5

20

6

19

7

18

8

17

A8

1

A2

3

22
21

A9

1

23

A10

A3

10

9
16

6

11

A4

5

3

10
15

11

14

13

12

A5

A7
A6

8

9

(a) Circular-arc model CM7(b) Circular-arc graph G

Figure 1: Circular-arc model CM and graph G
Table 1: Details of an circular-arc model CM.
i
ai
bi

3.

1
20
1

2
24
4

3
5
6

4
2
7

5
3
10

6
8
13

7
9
14

8
12
17

9
15
18

10
16
21

11
11
22

12
19
23

Algorithm DHV

In this section, we present Algorithm DHV (Detour Hinge Vertex) for solving the detour
hinge vertex problem of circular-arc graphs. We describe some definitions and lemmas
which are useful for constructing the algorithm for the problem. The detail descriptions of
proofs are omitted due to limitations of space.
The set of the minimum number of arcs that covers the entire circumference C of CM isdefined as MCA. The shortest circumference circuit of circular-arc graph G corresponding to
MCA is denote as S CC. In an example of Figuer 1(a), we have MCA = {A1 , A2 , A5 , A6 , A11 }and S CC = ⟨1, 2, 5, 6, 11⟩. Note that there may be multiple shortest circumference circuits
in G. In Figuer 1(a), G has two shortest circumference circuits S CC1 = ⟨1, 2, 5, 6, 11⟩
and S CC2 = ⟨1, 2, 5, 7, 11⟩. The vertex commonly included in all shortest circumference
circuits is called a duplicate vertex. In this example, vertices 1, 2, 5, and 11 are duplicate
vertices. Let U be a hinge vertex set of G. We define U1 = {u | u ∈ U, u ∈ S CC} andU2 = {u | u ∈ U, u < S CC}. For S CC1 = ⟨1, 2, 5, 6, 11⟩, we have U1 = {1, 2, 5, 11} andU2 = ∅.
Lemma 1. For a circular-arc graph G = (V, E), v < V is not adjacent to discontinuity two
vertices in S CC of G.
Lemma 2. For a circular-arc graph G = (V, E), v ∈ V is adjacent to at most three vertices
in S CC of G.
Lemma 3. For a hinge vertex u < S CC in a circular-arc graph G, we have det(u) ≤ 4.

31

Lemma 4. For a hinge vertex uinS CC in a circular-arc graph G, det(u) ≤ 4 if u is not
duplicate vertex.
Lemma 5. The det(u) of hinve vertex u ∈ U1 can be computed in O(n2 ) time.
Lemma 6. The det(u) of hinve vertex u ∈ U2 can be computed in O(n2 ) time.
Algorithm DHV can solve the detourl hinge vertex problem of a circular-arcs graph G.
Algorithm 1: Algorithm DHV
Input: All terminal points ai , bi for n circular-arcs in M.
Output: The detourl hinge vertex of a circular-arcs graph G.
(Step 1) /* Computation of hinge vertex set U */Construct an extended circular-arc model ECM from M ;Comput hinge vertex set U of G by applying Honma et.al’s algorithm [6] ;(Step 2) /* Compute MCC, U1 , U2 */Compute the shortest circumference circuit MCC ;Set U1 := {u | u ∈ U, u ∈ MCC} ;Set U2 := {u | u ∈ U, u < MCC} ;(Step 3) /* Compute detour degree det(u) */for each u ∈ U1 do compute det(u) ;for each u ∈ U2 do compute det(u) ;(Step 4)/* Maximum detour hinge vertex of G */Maximum detour hinge vertex is k such that det(k) = max{det(u) | u ∈ U}.

4.

Conclusion

In this paper, we proposed an algorithm for solving the detour hinge vertex problem ofcircular-arc graphs in O(n2 ) time. Reducing the complexity of the algorithm and extending
the results to other graphs will be considered in future research.

References[1] J.M. Chang, C.C. Hsu, Y.L. Wang, and T.Y. Ho, “Finding the set of all hinge vertices for strongly chordal
graphs in linear time,” Inf. Sci., vol.99, no.3-4, pp.173–182, 1997.
[2] H. Honma and S. Masuyama, “A parallel algorithm for finding all hinge vertices of an interval graph,”
IEICE Trans. Inf. & Syst., vol.E84-D, no.3, pp.419–423, 2001.
[3] H. Honma, Y. Nakajima, Y. Igarashi, and S. Masuyama, “Algorithm for finding maximum detour hingevertices of interval graphs,” IEICE Trans. Fundamentals, vol.E97-A, no.6, pp.1365–1369, 2014.
[4] M. Pal, “Intersection graphs: An introduction,” Annals of Pure and Applied Mathematics, vol.4, no.1,
pp.43–91, 2013.
[5] M.C. Golumbic, Algorithmic Graph Theory and Perfect Graphs, Academic Press, New York, 1980.
[6] H. Honma and S. Masuyama, “An optimal parallel algorithm for finding all hinge vertices of a circular-arc
graph,” IEICE Trans. Fundamentals, vol.E91-A, no.1, pp.383–391, 2008.
[7] F.R. Hsu, K. Shan, H.S. Chao, and R.C. Lee, “Some optimal parallel algorithms on interval and circular-arc
graphs,” J. Inf. Sci. Eng., vol.21, pp.627–642, 2005.

32

JSST 2022

Development of a System to Predict Child Abuse in Newborns ∗Yuya Sato1 , Yoko Nakajima2 , Hirotoshi Honma2,* , Yasuko Makita3 , Kazuyo Matsuura3
1 Advanced Course of Electronic and Information System Engineering, National Institute ofTechnology Kushiro College
2 Department of Creative Engineering, National Institute of Technology Kushiro College
3 School of Nursing, Sapporo City University∗ honma@kushiro-ct.ac.jp

Abstract. The purpose of this study is to build a system that determines the risk of newborn
child abuse based on parents’ and newborns’ data (face data/observed data) collected by the
NICU staﬀ.
Keywords: predict child abuse in newborns, natural language processing, machine learning

1.

Introduction

The number cases involving child abuse in Japan have been increasing in recent years andis already considered a social problem. The number of children counselled in child consultation centers in Japan rose from 1,101 in 1990 to 122,575 in 2016 [1]. As the numberof child abuse cases continues to grow, abuse-related deaths are also becoming increasingly
common.
Previous studies have indicated that many cases of child abuse include newborns withlow birth weight and congenital diseases. Low birth weight of babies may be associated withsocial problems such as young women aiming to lose weight, increasing cases of late marriage/advanced maternal age, multiple pregnancies, negligence in prenatal health check-up,
as well as choosing an early delivery date in high-risk pregnancy cases, thanks to advancesin the management of newborns [2]. Moreover, newborns with congenital diseases accountfor 4 to 5 ％ of all births [3].
Premature babies, low birth weight babies, babies with congenital diseases, and those
who require resuscitation at delivery and are admitted to the NICU(Neonatal Intensive Care
Unit) for treatment. Many parents of newborns admitted to the NICU are highly likely
to commit child abuse due to a series of problems such as a higher economic burden ofparenting, the mother’s health problems, socio-economic problems, and unconducive homeenvironments. As the newborns admitted to NICU are in life-threatening situations, theparents are more susceptible to anxiety and depression. Moreover, if the parents are not∗

This work was partially supported by JSPS KAKENHI Grant Number 19K11834, 20K11093, and22k10924.

33

given appropriate support from the time their newborn is admitted to the NICU, their stresslevels remain or even increase, making them even more likely to commit child abuse after
the infants are released from hospital care. Hence, it is important that parents with newborns
in the NICU receive proper support and care to prevent them from harming their children.
Skilled NICU doctors and nurses can predict individual cases of newborn child abuse
to an extent. However, since they rely on only their instincts, it is possible to misjudge
the situation and end up overlooking the risk of abuse. Hence, the purpose of this study
is to build a system that determines the risk of newborn child abuse based on parents’ andnewborns’ data (face data/observed data) collected by the NICU staﬀ.

2.

Factors of Abuse

The main factors of newborn abuse are as follows [4]:1. Risk factors on the part of parents (unwanted pregnancy/childbirth, young pregnancy/childbirth,parents with mental illness, etc.)
2. Risk factors on the part of the child (premature birth, disabilities, and other factorsthat make parenting diﬃcult)
3. Risk factors related to child rearing (financial distress and social isolation of the family, single-parent families including unmarried couples, families with unstable relationships, families without support from those around them, families with frequent
relocations, husband-wife fights, violence from the spouse, etc.)
4. Other risk factors (not notifying pregnancy or birth registration, delay in application
for child healthcare, not undergoing maternity health check-up, negligence in maintenance/promotion of health of unborn child and mother, infrequent visits for maternity
check-up, lack of communication with healthcare professionals, etc.)

3.

The Flow of the Study

Face data represents the parents’ attributes (e.g., the number of children, their economic
stability/status, chronic diseases, history of abuse) quantified in the simplest way possible.
Observation data comprises records maintained by nurses on the parents’ behavior and attitude, from the moment the child is born until their admission to the NICU, as well as theparents’ communication with doctors and other nurses.
Face data is composed of the following items:
Mother’s information “Marital status”, “Physical illness”, “Mental illness”, “Intellectual
disability”, “History of abuse”, “History of domestic violence”, “Pregnancy history”,“Birth history”, “Childcare history”, “Economic status”
Information about pregnancy “Frequency of health check-up”, “Maternity record book”

34

Information about delivery “Delivery method”, “Delivery location”, “Illnesses”
1. Processing observed data
Before extracting information from the observed data, a cleaning process was initiated to improve learning accuracy. First, to counter the problem of inconsistency,
the terminology was standardized using a dictionary. This process aggregates multiple synonymous words into a single term. Unnecessary stop words were also removed. All nouns and verbs were extracted from the preprocessed observed datausing a morphological analyzer, and the respective occurrence frequency and TF-IDF(Term Frequency-Inverse Document Frequency) value were calculated.
2. Correct answer label
First, a document was prepared to state whether the parents analyzed the committed
abuse after the child was released from the hospital (or whether the nurse felt a riskof child abuse). Then, the document was compared with the face data sets and observed data as correct answer labels. With this, this system builds a model that canaccurately predict the correct answer label from each piece of face data and observed
data through supervised learning.
3. Analysis method
Because the features obtained from face and observed data (volume of vocabulary)
are huge, it is diﬃcult to use these date as training data for machine learning. Therefore, a dimension reduction process using a feature selection method and principal
component analysis was carried out. Then, we built three models: (1) using only theface data, (2) using only the observed data, and (3) using both sets of data. The respective prediction accuracy of each model was compared. Multipul classifiers wereused to compare the data models and verify the risk classification accuracy. To prevent overlearning in the calculation of accuracy, cross-validation was performed. The
assessment of this study focuses on accuracy and, especially, on recall. This is because the system is designed to detect the risk of child abuse that is, since it concerns
the life of newborn babies, absolutely nothing can be overlooked.

4.

ResultWe carried out an experiment to classify whether the nurse felt at risk or not. We used thescikit-learn module as a machine learning tool for classification. We adopted SVM, Logistic
Regression, Random Forest, and AdaBoost Classifier as classification learning models. 5
cross-validation performed to each classification results. Table 1 shows each classification
results.

5.

Conclusion

In this study, the risk of child abuse was determined using data from 166 cases.AdaBoostClassifier shows the highest value in Accuracy, Recall, Precision, F-value performance met-

35

Method
SVM
Logistic Regression
Random Forest
AdaBoost Classifier

Table 1: Classify result
Accuracy Recall Precision
0.778
0.500
0.390
0.778
0.500
0.390
0.881
0.732
0.934
0.995
0.988
0.997

F-value
0.437
0.437
0.779
0.992

rics. AdaBoost is a very popular boosting technique that aims at combining multiple weakclassifiers to build one strong classifier. Since the classification result of this system is related to life rescue, it is very important that the recall value is high.
The Model is still unreliable due to small number of data. However, we believe it can bea useful tool to improve the quality of medical care. If the algorithm used in this study canbe incorporated into electronic medical record systems, it would be possible to quantify andregister the impressions and risk gauged by medical staﬀ in every childbirth case, and canalso be shared much more eﬃciently. Also, if the staﬀ realizes that their notes help directly
identify the risk of child abuse, it may make them more sensitive to this problem and write
their reports more accurately.

References[1] 厚生労働省, “児童相談所での児童虐待相談対応件数とその推移,” 厚生労働省第 20回児童虐待防止対策協議会, 2017.
[2] 中村肇, “超低出生体重児の育ちを見守る,” 学術の動向, vol.15, no.4, pp.15–21, 2010.
[3] 仁志田博司, 新生児学入門, 医学書院, 東京都, 2018.
[4] 厚生労働省雇用均等・児童家庭局総務課, “子ども虐待対応の手引き,” 2013.

36

JSST 2022

A System to Generate Kanji Reading and Writing
Questions for Learners of the Japanese Language∗
Toshiki Aisaka1 , Yoko Nakajima1 , Hirotoshi Honma1,* , Tomoko Watanabe2
1 Department of Creative Engineering, National Institute of Technology Kushiro College
2 Graduate School of Humanities and Social Sciences, Hiroshima University∗ honma@kushiro-ct.ac.jp

Abstract. In this study, we developed a question generation support system for kanji reading
and writing for Japanese learners. This system adopted a multiple-choice question method.
In addition to helping instructors develop reading comprehension tests, this system can beused to support autonomous learning for students of Japanese.
Keywords: natural language processing, generation of kanji reading and writing exercises

1.

Introduction

Many Japanese-language education institutions inside and outside Japan provide reading
comprehension instruction to students. Computers have been used in Japanese-language
reading comprehension instruction since the 1980s, primarily to provide learning materials
and select teaching materials from the Internet [1, 2]. However, the authors of this paper
were unable to find research focused on the use of computers in reading comprehension
assessment.
Reading comprehension is commonly assessed through specific tests designed to verify
students’ ability to understand given texts. Instructors spend a great deal of time and eﬀort
writing the questions of these tests, either by reusing questions used during class or by
writing questions based on new but similar texts.
We are currently building a system that automatically generates reading comprehension
questions for Japanese texts in order to support students (especially international students)
who are learning Japanese. We have already developed systems that automatically generate
supplementary noun/verb phrase questions, conversational sentence rearrangement questions, and supplementary questions on frequently appearing words.
In this study, we developed a question generation support system for kanji reading andwriting for Japanese learners. In addition to helping instructors develop reading comprehension tests, this system can be used to support autonomous learning for students of Japanese.
∗

This work was partially supported by JSPS KAKENHI Grant Number 19K11834, 19H01270, and20K11093 and Nagaoka University of Technology (NUT) grant for collaborative research with National Institute of Technology (NIT).

37

2.

Kanji Reading Questions

Kanji reading questions are answered by choosing the correct reading of a kanji in a specified part of a sentence. Learners choose from the given answer options. Each answer option
includes four possible answers: one correct answer and three wrong answers called “distractors.” Learners can specify the diﬃculty level of the question. The system automaticallygenerates questions with varying levels of diﬃculty that align with the language-ability levels N1 to N4 in the Japanese-Language Proficiency Test (JLPT).
We describe the process of automatic generation of kanji reading questions.
The question creator decides which kanji to use as an answer to the question. The targetof the reading questions and the kanji to be answered are determined from among the kanjiacquired in the lower grades of Japanese elementary schools for learners at the N4 level, inthe upper grades for learners at the N3 level, in Japanese junior high schools for learners at
the N2 level, and in high schools and above for learners at the advanced N1 level. However,the parts of speech of kanji are limited to nouns, adjectives, and verbs. Naturally, many of
the kanji learned in the lower grades of elementary school are basic and appear with highfrequency in Japanese texts.
Texts containing the selected kanji are automatically extracted from various corpora(newspapers, novels, fairy tales, nursery rhymes). The extracted texts are automaticallysummarized so that they are not longer than seven sentences. The summarization mechanism
uses MeCab and CaboCha tools, and unnecessary modifications are removed based on SVO
syntax. If the summarized text contains one or more kanji of a higher diﬃculty level thanthat assumed by the question creator, the text is discarded.

2.1.

Generating Distractors for Kanji Reading Questions

Each answer option includes four possible answers: one correct answer and three wrong
answers called distractors. Many studies have been conducted on the automatic creation ofanswer options [3, 4]. We describe our method for generating distractors.
[R1] The first method is to include among the distractors a diﬀerent reading from thecorrect one. For example, for a kanji whose correct answer is the on-yomi reading, the kunyomi reading is included in the distractors. Similarly, when the correct answer is a kanji
with a kun-yomi reading, its on-yomi reading is included in the distractors.
[R2] As a second method, the readings of kanji that are synonyms of the correct answer
or have an opposite meaning to the correct answer are included among the distractors.
[R3] For questions in which students are asked to read a kanji with a jyukugo-kun reading that is neither an on-yomi reading nor a kun-yomi reading, the reading of an associative
word of the kanji is included in the distractors.

3.

Kanji Writing Questions

Kanji writing questions are questions to assign the appropriate kanji to specified hiragana
parts in a text. The structure of the answer options is the same as in the kanji reading
questions.

38

The question creator decides which kanji to ask the learners to write. Text containing
the target kanji is automatically extracted from various corpora and summarized so that thelength is no more than seven sentences.

3.1.

Generating Distractors for Kanji Writing Questions

As in the case of kanji reading questions, there are four answer options for kanji writing
questions, which consist of one correct answer and three distractors.
[W1] The first method is to include among the distractors other kanji that have the same
reading as the correct answer. This process is implemented using a tan-kanji dictionary.
[W2] In the second method, we include among the distractors kanji that are similar tothe correct answer or that easily come to mind based on the correct answer. This process isrealized using the kanji synonym dictionary and associative kanji dictionary.
[W3] The third method includes homophones of the correct kanji among the distractors.
There is a large number of homophones in kanji, especially in jyukugo (idiom). This process
is implemented using a homophone dictionary.
[W4] In the fourth method, kanji similar in shape to the correct answer are included
among the distractors. Since kanji belonging to the same radical class are similar in shape,
other kanji with the same radical as the correct answer can be used as distractors.
The questioner inputs the diﬃculty level of the question and the correct kanji into thesystem. The system automatically generates multiple kanji reading and writing questions
according to the given diﬃculty level, and multiple answer options corresponding to thosequestions.
The questioner brushes up to a more sophisticated question by selecting each generated
question and answer option, and modifying it by editing. Completed questions are stockedin the question database according to the diﬃculty level.
Learners can expect to improve their Japanese proficiency by inputting the appropriatediﬃculty level and question type (reading / writing question) into the system and solving thequestions given.

4.

Verification Experiment

We verified the system’s validity and accuracy by monitoring three Japanese (Japanese A,
Japanese B, Japanese C) and five international students. Among the international students,
one had Japanese proficiency at the JLPT N1 level, three at the N2 level (N2 level A, N2
level B, N2 lebel C), and one at the N3 level. Table 1 lists the results of tests conductedusing the 20 supplementary conjunction questions generated by the system. The first column of the table shows the problem tester axis. The second column of the table shows thequestion number, and the third column shows how many of the 20 questions were answered
correctly. Correct answers are indicated by a circle(⃝), and incorrect answers are indicatedby a cross(×).
All the Japanese participants achieved perfect scores, while the ratio of correct answers
among the international students varied based on their individual Japanese-language ability. The test participants rated their satisfaction with the test as generally high. Going39

forward, we plan to generate more diverse questions for the tests, as well as to implementa mechanism that automatically generates questions of diﬃculty-levels compatible with theJapanese-language abilities of the test-takers.
Table 1: Results of verification experiment
TesterJapanese A
Japanese B
Japanese C
N1 levelN2 level A
N2 level B
N2 level C
N3 level

5.

1
⃝
⃝
⃝
⃝
⃝
⃝
⃝
×

2
⃝
⃝
⃝
⃝
⃝
×
⃝
⃝

3
⃝
⃝
⃝
⃝
×
⃝
×
×

4
⃝
⃝
⃝
×
×
⃝
×
×

5
⃝
⃝
⃝
⃝
⃝
⃝
×
⃝

6
⃝
⃝
⃝
⃝
⃝
×
⃝
×

7
⃝
⃝
⃝
×
⃝
⃝
×
⃝

8
⃝
⃝
⃝
⃝
×
×
×
×

Question number
9 10 11 12 13⃝ ⃝ ⃝ ⃝ ⃝⃝ ⃝ ⃝ ⃝ ⃝⃝ ⃝ ⃝ ⃝ ⃝× ⃝ × ⃝ ⃝× ⃝ ⃝ × ⃝× ⃝ × ⃝ ⃝⃝ ⃝ ⃝ × ⃝× ⃝ × ⃝ ×

14
⃝
⃝
⃝
⃝
⃝
⃝
⃝
⃝

15
⃝
⃝
⃝
⃝
⃝
⃝
⃝
⃝

16
⃝
⃝
⃝
⃝
⃝
⃝
⃝
⃝

17
⃝
⃝
⃝
⃝
×
⃝
⃝
×

18
⃝
⃝
⃝
⃝
×
⃝
×
⃝

19
⃝
⃝
⃝
×
×
×
⃝
×

20
⃝
⃝
⃝
⃝
⃝
×
⃝
×

Correct
20
20
20
15
12
13
13
9

Conclusion

In this study, we developed a system to support the creation of kanji reading and kanji
writing questions. The quality and accuracy of the questions generated by the system werelargely evaluated favorably in the verification experiments.
When the question creator inputs the diﬃculty level of the question and the correct
kanji, the system automatically generates kanji reading and writing questions according to
the given diﬃculty level, as well as multiple answer options corresponding to the questions.
The question creator sifts through the generated questions and answer options and edits
them as necessary to polish and refine the questions. The completed questions are stored in
the question database according to their diﬃculty levels. This system significantly reduces
the amount of work required and the burden on the question creator and can dramaticallyincrease the productivity of question writing.
Learners can input the level of diﬃculty and the type of question (reading or writing
questions) according to their own proficiency level into the system, and solve the questions
to improve their Japanese-language ability.
Our future challenges include the automatic generation of more appropriate questions
that match the learners’ level of proficiency, oﬀering questions from fields and genres ofinterest to the learner, and implementing the system as a web application.

References[1] A. Fukusaka, “Practice and evaluation of a problem producing exercise by using an automatic question
generation system,” JSiSE Research Report, vol.32, no.5, 2018.
[2] M. Furudate, “A proposal and evaluation on a method of automatic construction of knowledge base for
automatic generation of exam questions,” IPSJ SIG Technical Report, vol.2015-CE-128, no.14, 2015.
[3] T. Goto, “Automatic generation system of multiple-choice cloze questions and its evaluation,” Knowledge
Management & ELearning: An International Journal, vol.2, no.3, pp.210–224, 2010.
[4] T. Iwata, “Recommendation for english multiple-choice cloze questions based on expected test scores,”International Journal of Knowledge-Based and Intelligent Engineering Systems, vol.15, no.1, pp.15–24,
2011.

40

JSST 2022

脱塩基部位からなるクラスターDNA 損傷の分子動力学シミュレーション
Molecular Dynamics Simulation of Clustered DNA
Damage Composed of Apurinic/Apyrimidinic Sites
寺川 和志 1,2*, 藤原 進 1, 水口 朋子 1, 米谷 佳晃 2, 鹿園 直哉 2, 赤松 憲 2,
中村 浩章 3,4Kazushi Terakawa1,2*, Susumu Fujiwara1, Tomoko Mizuguchi1, Yoshiteru Yonetani2,
Naoya Shikazono2, Ken Akamatu2, Hiroaki Nakamura3,4

1

京都工芸繊維大学, 2 量子科学技術研究開発機構, 3 名古屋大学, 4 核融合科学研究所
Kyoto Inst. of Tech., 2QST, 3Nagoya Univ., 4NIFS

1*

m1672015@edu.kit.ac.jp

Abstract. DNA は酸化、放射線照射などに弱く、損傷する。脱塩基部位(AP 部位)、一本鎖切断（SSB）、二本鎖切断(DSB）といった DNA 損傷が、局所的に発生しクラスター化すると、細胞中の酵素による修復が困難となり、細胞死やがん化の原因になる。
本研究では、SSB や DSB につながると DNA に大きな影響を及ぼす可能性のある AP
部位に注目した。
AP 部位の損傷 DNA について分子動力学シミュレーションを行い、クラスター損傷と単独損傷の場合について、DNA の構造変化を解析した。また、分子構造の解析から損傷の生体影響について議論した。
Keywords: molecular dynamics simulation, apurinic/apyrimidinic (AP) site, clustered DNA
damage1. Introduction
Deoxyribo Nucleic Acid (DNA)はデオキシリボースとリン酸からなる主鎖と塩基から構成される多くの生物の遺伝情報の継承と発現を担う物質である。２本のポリヌクレオチド鎖は、相補的な塩基対(A-T 塩基対、G-C 塩基対)間の水素結合により結合

41

JSST 2022
し、右巻き二重らせん構造をとっている。
DNA はアルキル化や酸化、放射線照射などに弱く、その影響で損傷することがある。特に 2011 年の福島第一原子力発電所事故で、放射線物質汚染が問題となったが、放射線は細胞に当たると、細胞の中にある遺伝子の本体である DNA に傷をつけることがある。
DNA 損傷には、DNA のグアニン塩基が酸化した 8-オキソグアニン(8oxoG)、
塩基が脱離した脱塩基部位(apurinic/apyrimidinic site, AP 部位)、DNA 二重らせんの主鎖の片方が切れる一本鎖切断(single strand break, SSB)、または両方が切れる二本鎖切断(double strand break, DSB)などがある。これらの損傷が局所的に集まることもあり、クラスター損傷と呼ばれている。単独に生じた損傷は、通常、体内の酵素により修復されるが、傷が多かったり、クラスター化すると酵素による修復が不十分になり、細胞死やがんに繋がることがある。
放射線損傷による DNA 分子構造への影響は、分子動力学(MD)シミュレーションにより研究されてきた[1,2]。AP 部位が相補鎖内に２つあるクラスター損傷については調べられてきたが、同一鎖にある場合については研究されていない。両方の AP 部位
が体内の酵素によって別々に修復されると考えられ、重要視されてこなかったが、その修復や生物影響はよく分かっていない。本研究では、同一鎖上の AP クラスターについて MD シミュレーションを行った。MD 計算は DNA の構造情報を分子レベルで導くことができる。MD 計算からクラスター損傷によって DNA の構造がどのように変化するのかを明らかにした。また、同一鎖の AP クラスターと相補鎖の AP クラス
ターの違いも明らかにする。

2. Models & Simulation
以下のように未損傷 DNA と損傷 DNA に対して MD シミュレーションを行った。
シミュレーションソフトは、分子設計、MD 計算、計算結果の解析などができる、
AmberTools20 を用いた[3]。
2.1 未損傷 DNA
AmberTools20 の Nucleic Acid Builder モジュールを用いて未損傷の DNA を作成した。
塩基配列は鹿園らが実験で使用したものの一部を用いた[4]。次に、AmberTools20 のxLEAP モジュールを用いて、電荷を中和するために対イオン (Na+)を入れ、DNA の周りを溶媒(水)で満たした。DNA の力場は AmberTools20 の bsc1 力場を、水の力場は
TIP3P 力場を用いた。切頂八面体の周期境界条件(両端の距離は約 90Å)を設定した。
続いて、溶媒と対イオンのエネルギーを最小化した後、系全体のエネルギーを最小化した。得られた配置を初期配置として、周期境界ボックスの体積を固定して MD 計算を開始した。20ps(10,000 ステップ, 1 ステップ = 2 fs)かけて、系の温度を 0 K から 310
K まで昇温した後、温度を 310 K に維持し、一定圧力(1 atm)で 9ns の平衡化を行った。
2.2 損傷 DNA
このようにして平衡化した未損傷 DNA(Fig.1 (a)) に、AP 部位を１か所または２か所挿入した(Fig.1 (b))。AP 部位の原子電荷は、Gaussian 電子状態計算から求めたものを用いた。上記未損傷 DNA の場合と同様にして、xLEAP モジュールを用いて、対イオン (Na+)、溶媒(水)を加え、エネルギー最小化を行った。周期境界ボックスの体積を固定して MD 計算を開始した。系の温度を 0 K から 310 K まで昇温した後、温度を310 K に維持し、一定圧力(1 atm)で 50ns の MD 計算を行った。上記のプロセスを Fig.2に示す 7 種類の損傷 DNA に対して行った。損傷 DNA の構造変化を調べるため、損傷 DNA の根平均二乗偏差(root mean square deviation, RMSD)と C1’- C1’距離を解析した。比較のため未損傷 DNA に対しても同様の解析を行った。

42

JSST 2022
(a)

(b)

Fig.1 Snapshots of (a) an undamaged
equilibrated DNA (“undamaged” in Fig.2)and (b) a damaged DNA (“seq-5” inFig.2).

Fig.2 Models of clustered DNA damage composed of AP sites: an undamaged DNA, damaged
DNA with one AP site (seq0) and damaged DNAs
with two AP sites in the same strand (seq±1, 3, 5).
X represents an AP site.

Fig.3 Temporal evolution of RMSD of damaged DNAs composed of AP sites and an undamaged
DNA.

3. Results
3.1 Global structural change
AP 部位による構造の変化を解析するため、RMSD を計算した。RMSD は以下の式で表され、任意の構造｛𝒓𝑖 ｝と参照構造｛𝒓0𝑖 ｝との差を記述する。これは構造がどの程度変化したかの指標となり、RMSD の値が大きいと、構造の変化が大きいといえる。
𝑖は計算した原子の番号を表し、𝑁は計算した原子の個数を表す。
𝑁

RMSD = √∑𝑖=1(𝒓𝑖 − 𝒓0𝑖 )2𝑁

DNA 主鎖（P, O, C を対象）についての RMSD を 2ps ごとに計算した。未損傷 DNA
及び損傷 DNA の時間に対する RMSD を Fig.3 に示す。Fig.3 によると、未損傷 DNAの RMSD(undamaged)より、損傷 DNA の RMSD(seq±1, 3, 5)の方が大きい傾向にあると判断できる。全時間の RMSD の平均値を比較すると、seq0 < undamaged < seq-3 < seq+1≈ seq-1 ≈ seq+3 < seq+5 < seq-5 となった。単独損傷(seq0)の場合は、未損傷(undamaged)の場合より若干小さいが、それほど変わらなかった。一方、複数損傷(seq-5)の場合では、主鎖の曲がりが他の場合より大きい。これが RMSD が大きい原因と考えられる。

43

JSST 2022

Fig.4 Temporal evolution of C1’-C1’ distance of damaged DNAs composed of AP sites and an
undamaged DNA.
複数の損傷が入ると未損傷の場合と比べて DNA の構造が全体的に変化しているといえる。
3.2 Local structural change
AP 部位による局所的な構造変化を解析するため、損傷（脱塩基部位）とその相補塩
基対のデオキシリボースの C1’原子間の距離を計算した。
各 DNA の C1’-C1’距離を Fig.4に示す。未損傷 DNA では、C1’-C1’距離が、10.69±0.17Åであった。損傷 DNA では、
C1’-C1’距離が増減しており、構造が変化していることがわかる。seq-1 では、AP 部
位の向かい側のシトシン塩基が DNA 二重らせんの外に配置している構造が確認できた。seq+3 では AP 部位が、らせんの外に位置しており、塩基対が上下に 1 ペアずれたような、チミン－チミンという異常な塩基対の形成が見られた。seq0 では AP 部位
が、らせん外に位置したことにより、チミンが、アデニンとシトシンの 2 つの相手と多重に塩基対を形成していることを確認できた。DNA の局所的な構造についても、損傷を受けると大きく変化することが分かった。このような構造変化が、酵素による損傷修復にどのような影響を与えるのかについて、今後検証していく必要がある。

4. Summary
AP 損傷が同一鎖上に生成した場合においても、DNA の分子構造は大きな影響を受けることが、分子動力学シミュレーションにより分かってきた。このような分子構造の変化が、酵素による損傷の修復に影響を与える可能性について検証していく必要がある。今後は、より長い時間の MD シミュレーションを行って DNA の構造変化をさらに解析し、損傷 DNA の構造変化と酵素による修復の関連性を明らかにする。

References
[1]
[2]
[3]
[4]

E. Bignon, H. Gattuso, C. Morell, F. Dehez, A. G. Georgakilas, A. Monari, and E.
Dumont. Nucleic Acids Research, 44 (2016), 8588-8599.
H. Fujimoto, M. Pinak, T. Nemoto, P O’Neill, E. Kume, K. Saito, H. Maekawa, Journalof Computational Chemistry, 26 (2005), 788–798.
D.A. Case, K. Belfon, I.Y. Ben-Shalom et al. (2020), AMBER 2020, University of California, San Francisco.
N. Shikazono, P. O‘Neill, Mutation Research, 669 (2009), 162–168.

44

JSST 2022

CGAN を用いた自動コード推定
Using CGAN Automatic Code Estimation
森田 龍斗 1,* , 原 崇 2
Ryuto Morita1,* , Takashi Hara2
1 鹿児島工業高等専門学校 専攻科 電気情報システム工学専攻
2 鹿児島工業高等専門学校 情報工学科
1 Advanced Electronic and Information Systems Engineering,

National Institute of Technology, Kagoshima College
2 Department of Information Engineering,
National Institute of Technology, Kagoshima College∗ b16166@kagoshima.kosen-ac.jp

本論文のフルペーパーは，日本シミュレーション学会論文誌に掲載されています．

DOI: https://doi.org/10.11308/tjsst.15.42

45

JSST 2022

深度画像による形推定を利用した道路標識の識別
Identification of Road Signs
Using Shape Estimation from Depth Images
愛甲 遼 1*, 武田 和大 2, 原 崇 2
Haruka Aiko1*, Kazuhiro Takeda2, Takashi Hara2

1

鹿児島工業高等専門学校 専攻科 電気情報システム工学専攻
2
鹿児島工業高等専門学校 情報工学科
1

Advanced Electronic and Information Systems Engineering,
National Institute of Technology, Kagoshima College
2
Department of Information Engineering,
National Institute of Technology, Kagoshima College*b16126@kagoshima.kosen-ac.jp

Abstract. Identification of road signs is important for automated driving of automobiles. In
some cases, road signs cannot be correctly identified because some of them are hidden by obstacles such as roadside trees, or because their patterns are similar to those of road signs, resulting in misrecognition. Therefore, this study aims to improve identification accuracy by removing obstacles using depth images. In addition, by understanding the shape of the object,misrecognition of objects with similar patterns is reduced.
Keywords: Automatic Driving, Road Sign, Obstacle, Identification

1. 背景と目的
近年，AI 技術は目覚ましく発展してきており，人々の生活をより豊かにしようと様々な分野に応用されている．
その例として自動車の自動運転化が挙げられる．自動車の運転には高い集中力と運転技術が必要とされ，自動運転化が実現することで運転自体への負担が大きく軽減される．
自動運転化を進める上で，道路標識を識別することが重要である[1]．ただし，道路標識
は街路樹などにより１部が隠れている場合や，お店の看板などの模様が道路標識と似ているため，誤検出をする場合など，道路標識の識別を阻害する要因がいくつかある．
そこで，本研究では深度カメラを用いて，次のような阻害要因を除去する手法を提案する．ひとつは，道路標識とその前景にある物体を検知し，障害物を除去すること．もうひ49

JSST 2022
とつは，形状を把握することで，模様が似ている物体の誤認識を低減することである．それらを実現することによって，識別精度の向上を図る．なお，識別に用いるニューラルネ
ットワークは，画像処理というタスクにおいて高い精度が報告されている[2]，畳み込み
ニューラルネットワーク（CNN）とする．

2. 提案手法
まず，深度カメラを用いて画像と深度データを取得する．今回は深度カメラとしてIntel の RealSense[3]を用いる．
次に，取得した深度データを用いて，各深度（距離）で画像を分割し，分割した画像から道路標識以外の前景や背景となる部分を除去する．
その後，色や形などから標識をグループ分けし，グループに応じた CNN へと入力する．
システム概要図を図 1 に示す．

図１．システム概要図

3. 予備実験
提案手法を行う前に，現在の CNN の性能について確認を行った．
データセットは CNN で学習するための道路標識画像である．識別する道路標識を図２に示す 24 種類とした．
データセットは，まず，道路標識が写っている画像から道路標識の部分をトリミングしたものを各標識について 10 枚ずつ用意した．次に実際に運転中の車上から撮影した画像に張り付ける．この際，道路標識のサイズは張り付ける画像の 1/8 から 1/10 の大きさとし，ランダムに張り付ける．また，張り付ける画像の場所は実際の車上から撮影した際に標識が現れる可能性が高い画像の左上部分に，ランダムに配置した．また画像には，明るさの変更，回転，ぼかし，ノイズなどを付加した．また，最終的な画像の枚数は速度制限の標識が各 6000 枚，それ以外の画像は各 5500 枚とした．速度制限の標識は，ほかの標識と比べて識別の重要性が高いため，画像枚数を多くした．
作成したデータセットをトレーニング用データとテスト用データに 8:2 の割合で分割する．まず，CNN にトレーニング用データを入力画像とし，学習させる．次に，テスト用
データを入力画像とし，作成したニューラルネットワークが正確に道路標識を分類できているか評価する．このとき，画像サイズは 169×169 とし，ネットワークもそれに対応した形状とする．
このトレーニング用データを学習し，テスト用データで評価する作業を繰り返すこと

50

JSST 2022
で，精度の高いニューラルネットワークを作成する．
また，作成したデータセットとは別に，車上から撮影した道路標識が写っている画像を用意し，これを検証用データとする．

図２．識別する道路標識

図３．CNN モデル

4. 予備実験結果
今回使用したモデルは，畳み込み層 6 層，プーリング層 3 層，全結合層 2 層の計 11 層からなるネットワークである（図３参照）．なお，図の畳み込み層の 3×3×32 は畳み込みフィルタのサイズが 3×3 であることを示し，
チャネル数が 32 であることを示している．
畳み込み層とプーリング層を交互に構成することで畳み込み層が道路標識の特徴を抽出
し，プーリング層で平均化することにより，より多くの道路標識の特徴を認識できるようになっている．
テスト用データと検証用データを用いた時の正答率を表１に示す．なおこの値はそれ
ぞれのデータを各 10 回テストした結果である．表１より，テスト用データは 99.7％，検証用データは 85.4％の精度で道路標識を認識できた．
正解した画像の例を図４に，誤答した画像の例を図５に示す．誤答した画像の特徴は，道路標識の色と風景の色が類似している，道路標識の色が雨などで劣化している，入力画像自体がぶれているなどが原因で誤答したことがわかった．
また，自動運転の実現には入力画像から道路標識をリアルタイムで認識する必要がある．そこで入力画像から道路標識を認識する時間を測定した．認識時間を表２に示す．
表１．正答率

表２．認識時間データ

正答率[%]

枚数[枚]

認識時間[sec]

テスト用データ

99.70

100

8.6485

検証用データ85.40

10.0864

(a)

(b)

(a)

(b)

図４．正答した画像

図５．誤答した画像

51

JSST 2022

5. 検証用データの準備
本研究の効果を検証するために，道路標識が障害物で隠れていることを考慮した画像を作成する．この画像は予備実験で用いた検証用データの画像に処理を加え，道路標識の1 部分を隠すことで作成する．障害物は街路樹を想定しているため，道路標識を隠す部分
は左側，左上部分とする．また，人が道路標識を正しく認識できる範囲で隠すこととする．
障害物で隠れていることを考慮して作成した画像を，図６(a)～(c)に示す．また，標識と誤認識する画像の例を図６(d)に示す．

(a)

(b)

(c)

(d)

図６．検証用データ

6. 今後の予定
今回の予備実験の結果から，道路標識と背景との色の差が明確であれば，正確に識別する可能性が高いことが分かった．また，今回の検証に用いた画像は標識の全体が写っているものであるため，これらが障害物などにより１部が隠れている場合には，さらに正答率は下がる見込みである．これは，前景の影響によるものであるため，深度カメラを用いて道路標識を前景や背景と分離した後に，形推定して識別することで，識別精度を向上させる．
また，模様が道路標識と似ているため，誤認識してしまう場合にも，形推定にて除外することができる見込みなので，誤認識する場合を減らすことができる．
今後は，検証用データを用いて，提案手法の効果を測定する．

参考文献[1] 池田 浩輔，孕石 泰丈：低品質画像における領域畳み込みニューラルネットワークを用いた道路標識認識，山梨講演会講演論文集，2017 ．

[2] 福井 宏，山下 隆義，山内 悠嗣，藤吉 弘亘，村瀬 洋：Random Dropout と Ensemble
Inference Networks による歩行者検出と標識認識，情報処理学会論文誌，Vol.57，No.3，
pp.910-921，2016．

[3] Intel®RealSenseTMTechnology
https://www.intel.com.au/content/www/au/en/architecture-and-technology/realsenseoverview.html，(閲覧日：2022/01/26)

52

JSST 2022

電子ミラーに映される映像の高画質化
Enhancement of Images Reflected on
Electronic Mirrors.
東山 翔馬 1*, 古川 翔大 2, 原 崇 2, 武田 和大 2
Shoma Higashiyama1*, Shota Furukawa2, Takashi Hara2, Kazuhiro Takeda2
1

鹿児島工業高等専門学校 専攻科 電気情報システム工学専攻
2
鹿児島工業高等専門学校 情報工学科
1

Advanced Electronic and Information Systems Engineering,
National Institute of Technology, Kagoshima College
2
Department of Information Engineering,
National Institute of Technology, Kagoshima College*b16156@kagoshima.kosen-ac.jp

Abstract. Recently, with the progress of computerization of automobile mirrors, there are
increasing opportunities to see images taken at night. Normally, an image taken at nightshows a bright light source such as a street light, and another background or a building looks
dark in the image taken at night, resulting in poor visibility. Therefore, we propose a method
to brighten the nighttime image naturally by using the histogram created from the pixel values. This method makes it possible to make a dark image brighter and easier to see as awhole.
Keywords: Night Images, Electronic Mirror, Visibility

1.

背景と目的

近年、自動車のルームミラーやサイドミラーの電子化が進んでいる中で、夜間に撮影した映像を目にする機会が増えている。
通常、夜間に撮影した画像は、街灯などの光源部分が明るく映り、そのほかの背景や建物などは暗く映ってしまう。自動車の電子ルームミラーや電子サイドミラーの映像の場合も同様で、他車のヘッドライトによる白飛びや、暗い部分の黒つぶれが発生してしまい、車両周辺の視認性が悪くなる。暗い画像の明るさを調整する画像処理の手法はいくつかあるが、十分な明るさが得られない、色味が変わる、白飛びや黒つぶれの発生など、不自然になるデメリットがある。また、白飛びや黒つぶれを軽減する手法の一つとして、HDR があるが、HDR は、露光が異なる複数枚の画像を組み合わせて一枚の画像を作成する手法のため、動画の撮影には適していない。

53

JSST 2022
そこで、本研究では暗い場所や夜間に撮影した画像を、明るさを調整することによって見やすくする手法を提案する。提案手法を用いることで、画像内の白飛び・黒つぶれを軽減し、自然な色味の画像が得られる。それにより、視認性を向上させ、夜間運転時のストレスの軽減や事故を減少させることを目的とする。

2.

Contrast Pair 法

Contrast Pair 法[1]とは、トーンカーブによって画像の V 成分を変更し、画像のメリ
ハリを強調する処理手法の一つである。
まず、画像を RGB 画像から HSV 画像に変換する。
次に、画像の V 成分から、ヒストグラムを作成する。最初に、隣接するピクセル間のの画素値(V 成分)の差を取り、その差がしきい値よりも大きい場合は、ヒストグラム内の、２つのピクセルの画素値の間の範囲の高さを、それぞれ 1 ずつ上げる。これを画像
内の全画素について行い、ヒストグラムを完成させる。
最後に、完成したヒストグラムの累積ヒストグラムをトーンカーブとして、画像の V
成分を変更する。Contrast Pair 法の一連の流れを図 1 に示す。

図 1 Contrast Pair 法によるトーンカーブ作成

3.

提案手法

提案手法では、Contrast Pair 法を改良した手法によって、夜間画像の高画質化を図る。Contrast Pair 法では、画像のメリハリを強調することができるが、全体的に暗い
画像を明るくすることはできない。そこで、Contrast Pair 法で作成した累積ヒストグ
ラムに、画像の V 成分の累積ヒストグラムを合成することで、新たなカーブを作成する。そして、作成したカーブに基づいて画像の V 成分を変更する(改良版 Contrast Pair法)。この改良によって、画像のメリハリを強調しつつ、白飛びや黒つぶれを軽減することができる。改良版 Contrast Pair 法の一連の流れを図.2 に示す。
次に、提案手法の一連の流れについて、以下に示す。
最初に、RGB 画像を HSV 画像に変換し、その V 成分から、改良版 Contrast Pair 法
を用いてトーンカーブを作成する。
続いて、HSV 画像の V 成分に対して以下の式を用いて、ガンマ補正を行う。
さらに、改良版 Contrast Pair 法によって作成したカーブに基づき、HSV 画像の V
成分を変更する。

54

JSST 2022
これらの処理を施した後、HSV 画像から RGB 画像に変換する。変換した後に、以下の式を用いて、画像内のヘイズを除去しつつ、色彩を強調する[2]。
𝑅 ′ (𝑥, 𝑦) =𝑅(𝑥, 𝑦) − 𝑤 × min⁡(𝑅(𝑥, 𝑦), 𝐺(𝑥, 𝑦), 𝐵(𝑥, 𝑦))
1 − 𝑤 × min⁡(𝑅(𝑥, 𝑦), 𝐺(𝑥, 𝑦), 𝐵(𝑥, 𝑦))𝐺 ′ (𝑥, 𝑦) =𝐺(𝑥, 𝑦) − 𝑤 × min⁡(𝑅(𝑥, 𝑦), 𝐺(𝑥, 𝑦), 𝐵(𝑥, 𝑦))
1 − 𝑤 × min⁡(𝑅(𝑥, 𝑦), 𝐺(𝑥, 𝑦), 𝐵(𝑥, 𝑦))𝐵′ (𝑥, 𝑦) =𝐵(𝑥, 𝑦) − 𝑤 × min⁡(𝑅(𝑥, 𝑦), 𝐺(𝑥, 𝑦), 𝐵(𝑥, 𝑦))
1 − 𝑤 × min⁡(𝑅(𝑥, 𝑦), 𝐺(𝑥, 𝑦), 𝐵(𝑥, 𝑦))

これが提案手法の流れである。また、今回は𝛾 = 0.8、𝑤 = 0.7で設定してある。

図.2 改良版 Contrast Pair 法によるトーンカーブ作成

4.

実験

本研究では、撮影された映像を自動車の電子ミラーに映すことを想定して、走行中の自動車で撮影されたドライブレコーダーの映像から評価用の画像を作成し、その画像に対して提案手法を適用した。照明が多く、比較的明るい場所で撮影された画像と、照明が少なく、全体的に暗い場所で撮影された画像の例を図.3、図.4 に示す。

5.

実験結果

図.3、図.4 に対して提案手法を適用した結果を図.5、図.6 に示す。図.5 を見ると、図.3 では画像全体に白いもやがかかっていたのに対して、出力画像ではもやが消え、看板や車、歩道の色が鮮やかになっており、周囲の確認がしやすくなっていることがわかる。また、図.5 では図.3 と比べて、看板の白飛びが軽減されており、文字が見やすくなっている。図.6 を見ると、図.4 に比べて周囲の車や建物の色は鮮やかになっているが、画像全体が暗く、歩道部分の黒つぶれが軽減されていない。また、画像内のもやも残っており、視認性の向上は見られない。
それぞれの画像の平均輝度をまとめたものを表 1 に示す。表 1 を見ると、処理後の画像の平均輝度値は上がっていることがわかる。しかし、図.4 と図.6 を目視で比較すると、視認性の向上はあまりみられないことから、平均輝度値以外の評価指標が必要だと考えられる。

55

JSST 2022
表 1 それぞれの画像の平均輝度
図.3
図.4
処理前
92.9
65.8
処理後

図.5

図.6

111.0

82.8

図.3 明るい場所で撮影された画像

図.4 暗い場所で撮影された画像

図.5 明るい画像に処理を加えた結果

図.6 暗い画像に処理を加えた結果

6.

今後の課題

今回、夜間に撮影された画像に対し、白飛びや黒つぶれを軽減しつつ、画像の明るさを調整する手法を提案した。提案手法では、処理前の画像にある程度の明るさがある場
合は、画像の視認性を向上させることに成功した。しかし、光源が少なく、明るさがない場合は、周囲の確認に十分な視認性は得られなかった。
今後は、評価用画像の収集や、暗所で撮影された画像に対する手法の改良、画像の視認性について評価できる指標の導入に注力していく。

参考文献
[1]
STEFFI AGINO PRIYANKA，HSIAO-JUNG TUNG, YUAN KAI
WANG：CONTRAST ENHANCEMENT OF NIGHT IMAGES，Proceedings of the2016 International Conference on Machine Learning and Cybernetics, Jeju,
South Korea, 10-13 July
K.He, J.Sun and X.Tang, : Single Image Haze Removal Using Dark
Channel Prior, IEEE Transactions on Pattern Analysis and Machine Intelligence,
Vol.33, No.12, pp.2341-2353, Dec. 2011

[2]

56

JSST 2022

Numerical Experiments of Swelling Objects
Interacting with Newtonian Fluids
Niku Guinea1,* , Daisuke Toriu2 , Satoru Ushijima2
1 Department of Civil and Earth Resources Engineering, Graduate School of Engineering,

Kyoto University
2 Academic Center for Computing and Media Studies, Kyoto University∗ guinea.niku.25w@st.kyoto-u.ac.jp

Abstract.
Deformable swelling particles interacting with Newtonian fluid was computed using
the combination of direct-forcing immersed boundary method and mass-spring model. The
swelling of the object was simulated by changing the natural lengths of the spring models. In
addition, the solid-solid interaction is treated by utilizing the distinct element method. Twocases of interaction between multiple swelling objects and Newtonian fluid are proposed.
As a result, it was shown that the basic behaviors of the swelling-deformable objects arereasonably calculated with the present method.
Keywords: Fluid-solid interaction, Swelling, Deformable object, Immersed boundary method,
Mass-spring model

1.

Introduction

The interactions between Newtonian fluids and swelling bodies with deformable properties
are introduced in this study. The combination of direct-forcing immersed boundary method(DF/IB) [1] and mass-spring model (MSM) [2] is proposed for the computation of the interactions of fluid and deformable solid. In addition, the swelling of the object is proposed bychanging the natural lengths of the spring models. The solid-solid interaction is treated by
setting CDS (contact detection sphere) along the solid boundary to detect particle collisions
for DEM computation. The Lagrangian (mass) points used in the MSM-DF/IB computations
share the same points for the placement of CDS used in the particle collisions. Furthermore,as the swelling computation is utilizing the natural spring length used in MSM, the couplingof the three methods in one computation is fairly easy.

57

2.

Numerical method

The fluid is assumed to be incompressible, and the momentum equations are given by∂ui ∂(ui u j )
1 ∂p
∂2 ui
+
=−
+ ν 2 + fi + λi ,∂t
∂x j
ρ f ∂xi∂x j(1)

where u j is the velocity component in x j direction in two-dimensional Cartesian coordinates,
t is time, xi is the component of the Cartesian coordinate system, ρ f is density, ν is kinematic
viscosity, and p is pressure. In addition, ui is the velocity component, fi and λi are theexternal and fluid-solid interaction forces in xi direction.
For the details of the computation of MSM-DF/IB, the reader are referred to the detailed
explanation by Guinea et al. [3]. The swelling of objects is calculated by increasing thenatural spring length used in the mass-spring model. As the value of the natural spring
length increases, the restoring force acts on the mass points and the distances between themalso increases. The swelling scheme is applied to the computation of the swelling of asingle hydrogel particle and gives a good agreement when compared to the experiment as
introduced in Guinea et al. [4]. In the proposed method, it is assumed that the volume
increase of the object is equal to the absorbed fluid. Additionally, for simplicity, the densityof the object is assumed to be constant in this study.

3.

Applications

3.1.

Sedimentation of a single particle in fluid

A swelling particle with initial diameter D0 = 0.25 and density ρ s = 1.25 is set inside a
2 × 8 area. The particle is dropped from the initial position of (1, 6) and falls freely with thegravity acceleration of −980. The particle swells from t = 0.30 until the maximum diameter
Dmax = 0.30. The surrounding sides are treated as non-slip boundaries. The fluid density ρ f
and kinematic viscosity ν are 1.0 and 0.1 respectively. The spring constant and coeﬃcientof restitution of mass-spring model k s and e s are 4 × 106 and 0.5. The number of fluid cells
used in the computation is 200 × 800, the number of Lagrangian points Nl is 70, and ∆t is
5 × 10−3 The calculated particle velocity in y direction is observed and compared with thecase of sedimentation of non-swelling particle with D = 0.25 and D = 0.30.
The changes of the particle diameter and the time history of the particle velocities can beseen on Fig. 1. The computation shows reasonable results between the three computations.
On the non-swelling particles, the terminal velocity of the particle with smaller diameter(blue line) is smaller, causing the particle to reach the bottom side later than the particle withlarger diameter (green line). The red line represents the velocity of the swelling particle.
Prior to t = 0.3, the velocity of the swelling particle is equal to the blue line where theparticle is identical in size. Afterward, the particle starts to swell and the velocity graduallyreaches the green line. In addition, due to the increase of the particle velocity on the swelling
particle, it reaches the bottom side before the particle with D = 0.25 and after the particle
with D = 0.3.

58

(b) Time histories of particle velocity(a) Changes of particles diameter

Figure 1: Comparison of diameter and particle velocities between the three cases

3.2.

Multiple objects transported in two-sided lid-driven cavity flow

Multiple deformable solid particles with various sizes are arranged inside the 2D computational area filled with fluid as it can be seen on Figs. 2 (a) and 3 (a). The computational
area is filled with fluid with kinematic viscosity ν = 0.01 and fluid density ρ f = 1.0. The
surrounding walls are treated as non-slip boundaries. The top and bottom walls move in horizontal direction with velocities utop,1 = 1.0, utop,2 = 0.0, ubottom,1 = −1.0, and ubottom,2 = 0.0.
The Reynolds number Re is 100. 16 particles with the initial diameter D0 varied between
0.1 and 0.15 with density ρ s = 1.5 are set up on the computational area. The spring stiﬀness k s of the particles are set between 0.56 and 0.88, with the coeﬃcient of restitutione s = 1.0. For the DEM computation, the values of spring constant kD,n in normal directionis 1.0 × 107 . Tangential contact forces is assumed to be 0 in this computation. The coeﬃcient of restitution eD of DEM computation is 1.0, and the radius of the contact detection
spheres is 0.167D. The number of fluid cells used is 100 × 100, while the number of mass
(Lagrangian) points Nl is 21. In addition, the time increment ∆t is 1.0 × 10−4 . Two cases arecomputed in this study. In both cases, velocities on top and bottom walls are stopped aftert = 10. Non-swelling particles are introduced in the first case, while swelling particles areon the second case. The particles on the second case are set to swell linearly with the slopeof 3 × 10−3 until the maximum size Dmax of 1.5D0 .
Figures 2 and 3 show the computation results on the non swelling and swelling particles.
The black lines represent the line connecting the mass (Lagrangian) points, and red lines
represent the line connecting the outer side of the contact detection spheres. Additional lines
are drawn connecting the surface and the center of the object to help observe rotation on theparticles. Due to the movement of the top and bottom wall, deformations can be observed
especially closer to the moving top and bottom wall where velocity is higher and when
the solid-solid and solid-wall contacts occur. After the velocities are stopped on t = 10,
we can see the particles move only with their remaining speed and finally stopped. The
trajectories of the center points of the particles for 10 ≤ t ≤ 20 can be seen on Fig. 4. It canbe observed that the non-swelling particles are moving further than the swelling particles
after the velocities are stopped. In addition, after the velocity is stopped, the particles on
the second case are still swelling. As a result, the particles push each other, causing slightmovement as it can be observed on the discontinuous lines on Fig. 4 (b).

59

(a) t = 0(b) t = 5(c) t = 10(d) t = 15(e) t = 20

Figure 2: Computation of non-swelling particles(a) t = 0(b) t = 5(c) t = 10(d) t = 15(e) t = 20

Figure 3: Computation of swelling particles

0.8

0.8

0.6

0.6
x2

1

x2

1

0.4

0.4

0.2

0.2

00

0.2

0.4

0.6

0.8

00

1

x1

0.2

0.4

0.6

0.8

1

x1

(a) non-swelling particles(b) swelling particles

Figure 4: Trajectories of non-swelling and swelling particles for 10 ≤ t ≤ 20 (redpoints mark the position of the particles at t = 20)

References[1] M. Uhlmann : An immersed boundary method with direct forcing for the simulationof particulate flows, J. Comput. Phys., vol. 209, 448–476, (2005)[2] D. Terzopoulos, J. Platt, A. Barr, K. Fleischer : Elastically deformable models, ACM
Siggraph Computer Graphics, vol. 21, 205–214, (1987)[3] N. Guinea, D. Toriu, S. Ushijima : Computational method for interactions between
deformable objects and fluid flows using immersed boundary method and mass springmodel, Transaction of Japan Society for Simulation Technology, (2021)[4] N. Guinea, D. Toriu, S. Ushijima : Eulerian-Lagrangian Approach for Interactions
between Fluids and Multiple Deformable Swelling Objects using Mass-Spring Model,
The 25th Applied Mechanics Symposium, 2B07-12-01, (2022)

60

JSST 2022

ラグランジュ流体シミュレーションの直線型プラズマ実験解析へ向けた改良
Improvement of Lagrange fluid simulation
for application to a linear plasma experiment
峰 壮鴻 1, 星野一生 1
Morihiro Mine, Kazuo Hoshino1

慶應義塾大学 理工学研究科

Graduate School of Science and Technology, Keio University*

mine@ppl.appi.keio.ac.jp

Abstract.
Controlling heat and particle loads on the divertor plate is one of critical issues for nuclear
fusion reactor. Lagrange－Monte-Carlo(LG-MC) Scheme has been developed to analyze the
SOL plasma in three-dimension and evaluate heat and particle loads on the divertor plate. In
this study, to analyze a liner plasma device and validate the code, model for a particle recycling (production of neutral particle and its ionization) has been implemented into the
one-dimensional LG simulation. The improved code successfully reproduces the qualitative
characteristics of a particle recycling.
Keywords: nuclear fusion, divertor, SOL, Lagrange－Monte-Carlo(LG-MC), particle recycling

1. 研究目的・背景
核融合反応を起こすためには, 燃料となるプラズマを高温・高密度状態で長時間閉じ込めることが必要とされる.トカマク型核融合炉では, 磁場を制御することによってプラズマを閉じ込め, 核融合反応を実現する. しかし, 拡散等により一部の高温プ
ラズマはコア領域から漏れ出し, 周辺領域(SOL)中の磁力線に沿ってダイバータ板へ

61

JSST 2022と集中する. そのため, ダイバータ板は高い熱・粒子負荷にさらされ, 溶融・損耗が深刻な課題となっている. そこで, SOL 領域内で高温プラズマの輸送を制御し, ダイ
バータ板到達前に十分に低温化させることが必要となる. 実験に加え, 数値解析によりこのような制御手法の研究が進められてきたが, 従来の SOL 領域内の数値解析は 1
次元もしくは 2 次元の解析が主流であった. しかし, 核融合炉の設計研究が進展し,実際の壁形状を取り入れた 3 次元解析が必要とされている.
SOL プラズマの 3 次元解析に向けて, 粒子法による流体方程式の解法であるラグラ
ンジュ-モンテカルロ(LG-MC)法, 及び LG-MC 法を用いた電磁流体コードが開発された[1]. 理論的な解析解に基づいた数値的な検証が行われたが, 実際の実験解析への適用はされていない.
これまでに開発された LG-MC コードは基本部分のみであり, 実験解析に適用するにはコード改良が必要である. そこで, 本発表では, LG-MC コードの実験解析適用を目指したコード改良について報告する.

2. 方法
今回は, LG-MC シミュレーションの実験との比較に基づく妥当性検証に向けて, (1) 実験解析に向けた装置体系の導入, 粒子種の設定, (2) 物理モデルの追加を行った. まず,(1)については, 形状が単純で計測に柔軟性を持ち, 数値シミュレーションの検証が比較的容易な直線型プラズマ装置(例えば NAGDIS-Ⅱ[2])を想定し, モデル体型を調整した.
従来の LG-MC コードは水素プラズマを対象としていたが, 実験に合わせ He を扱えるよ
うにコードを拡張した. 次に(2)として, ダイバータ領域のプラズマ輸送で重要となるリ
サイクリング過程(プラズマが壁で再結合することで中性粒子となり, 装置内でイオン化することで再びプラズマとなる過程の繰り返し)を系に合わせた形で実装した.

3. 結果
軸方向の長さを 2 m とした 1 次元体系におけるリサイクリング効果の検証について述
べる. 𝑥 = 0 mから流入したプラズマは𝑥 = 2 mに位置する壁に向かい, 表面再結合により中性粒子となる. その後, プラズマ中で再びイオン化する, いわゆるリサイクリング過程が起こる. このリサイクリングモデルを従来コードへ組み込み, リサイクリングの有無によるプラズマ密度分布を比較した.
Figure 1 に LG シミュレーションにおけるリサイクリングによる軸方向の密度分布の違いを示す. リサイクリングを考慮しない場合(リサイクリング率 0 ％)とリサイクリン
グを考慮した場合(リサイクリング率 99 ％)のプラズマ密度の状態を比較すると, リサイ
クリングを考慮した場合は系全体において密度が増加しており, リサイクリングによる密度増加が再現できていることがわかる. 今回の解析条件では壁近傍(𝑥 = 2 m)のプラズ
マ温度が低いため, 中性粒子の発生源である壁近傍よりも, より上流側(0 m 付近)でイオン化が起こっており, 壁近傍より上流側でより大きく密度が増加していることが確認できる.
また一般的に, リサイクリングが活発になると, プラズマ密度の増加に加えて, イオン化によりプラズマのエネルギーが減少するため, プラズマ温度は低下することが実験やシミュレーションにより確認されている[3]. 改良した LG コードでも, Figure 2 に示すよ
うに, リサイクリングによる温度低下を再現することができている.

62

JSST 2022

Fig. 1 Density profile. Red and green symbols correspond to the case without and withrecycling, respectively.

Fig. 2 Temperature profile. Red and green symbols correspond to the case without and withrecycling, respectively.

4. 結論と今後の展望
本研究では,LGコードへ直線装置に合わせた形でリサイクリングモデルを導入

63

JSST 2022
し,1次元シミュレーションにおいて定性的に妥当な結果を得た.今後はMCコードも含めた3次元シミュレーションへ拡張し,定量的に実験データとの比較を進めていく
予定である.

謝辞
LGコードの改良にあたり, カラム研究所(英)の大澤瞭子博士には多くの助言を頂きました.ここに感謝いたします.

参考文献[1] R.Tatsumi, et al, Plasma Fusion. Res. 15 (2020) 1403003.
[2] N. Ohno, et al., J. Plasma Fusion Res. 78 (2002) 1093.
[3] A. Hatayama, et al., J. At. Energy Soc. Jpn. 50 (2008) 378

64

JSST 2022

ミャンマーの森林破壊防止と保全活動についての検討
Consideration of "Preventing deforestation and Forest
preservation activities in Myanmar"

ミョーティリトゥ, 大塚

晃

Myo Thiri Thu, Akira Otsuka

事業創造大学院大学

事業創造研究科

Graduate Institute for Entrepreneurial Studies, Niigata, 950-0916 Japan
jma11212065mt@jigyo.ac.jp , otsuka.akira@jigyo.ac.jp
Abstract. Deforestation become one of the most common problems around the world.
Deforestation occurs due to various reasons including farming, logging, mining, andurbanization. The loss of trees can cause flooding, climate change, soil erosion, etc.
Deforestation occurs in all regions including Myanmar. The population is increasing rapidly
and we need more space to survive more, to cultivate and consumption. Human have
responsibility to stop deforestation. In this paper, we will consider about the deforestation
effects and search a few solutions to solve this current problem in Myanmar which include
consideration on deforestation-free rubber plantation.
Keywords: Deforestation, Forests, Deforestation-free Rubber Plantations, Activities

1. はじめに
森林破壊は生物多様性、地球の生態系、地球温暖化、気候変動に悪影響を及ぼす
世界的な課題の一つである。人間が森林破壊を引き起こす原因となっている。FAO（国連食糧農業機関）によると、世界の森林減少の原因として、農業の拡大が 8 割近くを占め、その他、道路やダムなどのインフラ建設、鉱業活動、都市化などが原因となっている。ミャンマーは、多くの国民が農業に従事しているため、森林破壊の問題に直面している一つの国として知られている。天然ゴム栽培も森林破壊の原因の一つである。ミャンマーの天然ゴム産業は約 65 万ヘクタールの栽培地を有しており、世界最大のゴム生産地の一つである。ミャンマーの天然ゴム生産は年々増加
し、2017 年に世界の天然ゴム生産量の 10 位に達した。一方、違法伐採もミャンマー
において重要な課題となっている。2019∼2020 年に違法伐採は合計 1,2370 トンまで

65

JSST 2022
［1］

上回っている 。
本稿ではミャンマーの森林破壊の影響と森林減少防止のため、良質な自然林を切り開くのではなく、荒廃した土地に天然ゴム栽培などを検討する。さらに、違法伐採に関する環境保全林業省の対策について検討する。

図１森林減少と森林劣化の原因

出典：www.wwf.or.jp

図１は、世界各地の森林減少原因と減少状況を述べている。森林破壊に直面している地域別と森林減少の原因を「主な原因、二次的な原因、小さな原因、原因ではない」、4 つの色で示している。主な原因として、大規模農業開拓が見られる。

2. 世界の森林破壊状況

図２

国別の原生林の減少

66

出典：http://news.mongabay.com

JSST 2022
図２は、世界的な国別の原生林の減少の状況を最大 19 カ国で示している。世界の原生林面積は 1990 年以降、8,000 万ヘクタール以上減少している。図２より、2020
年の世界の原生林の減少は、ブラジル（原生林消失 170 万 ha）、コンゴ民主共和国（DRC）（49 万 ha）
、ボリビア（27 万 7 千 ha）、インドネシア（原生林消失 27 万 ha）
が上位を占めた。ミャンマーは４万 ha の消失で 17 位となっている。ミャンマーでは、2001 年から 2021 年の間に、シャン州での森林破壊が１万 4700ha と最も多く、タニンターは 487kha、カチンは 462kha と合わせて、国全体の損失の 56%を占める。
日本は 2010 年に 17 万２千 ha の森林があり、国土の 71％を占めている。2021 年には 18.6kha が失われた。
図３は、アジア 5 か国の 2002 年から 2020 年の、原生林面積の減少をヘクタール
で示している。2020 年にはラオスが一位となっている。ヨーロッパに比較して、アジアの減少傾向が高い[2]。

図３

アジア 5 箇国の原生林の減少

出典：http://news.mongabay.com

上述したように、世界的な森林破壊は激しくなっている。森林破壊の防止対策として、木を植えることが最も有効である。その他、紙の使用量を減らすこと、リサ
イクルすることも重要である。しかしながら、本稿ではミャンマーの森林破壊の減少を目指し、ミャンマーにおける最も適切な保護活動について検討していく。

3. ミャンマーの森林保全活動
近年、ミャンマーでは植林活動が活発に行われ、森林面積が増加してきた。2015年から 2020 年の間に、森林減少の割合は 1990 年代の年間 1,600 万ヘクタールから、
1,000 万ヘクタールに減少したと推定される[2]。

67JSST 2022
WWF（World Wildlife Fund）ジャパンと WWF（World Wildlife Fund）ミャンマーの協力で、2017 年に森林保全と持続可能な天然ゴムの生産を目指した新たなプロジ
ェクトを開始した。プロジェクトは、南北の保護区をつなぐタニンダーリ地域において、ミャンマー政府と小規模農家、また、カレン民族同盟とも協力してプロジェクトを実施した。守るべき森林と農園の開発適地を明確に区分けするための土地利用計画の作成を目指している。

出典：www.wwf.org.mm
また、環境や人権に配慮しながら、天然ゴム産業を振興するための新たな法律の制定も進めている。さらに、地域の小規模農家に対しては、天然ゴムの生産性と品質を向上させることを目的としたトレーニングの機会を提供する取り組みも行なっている。WWF は政府の農業技術指導員や小規模農家に対して、森林破壊を引き起こす
ことなく、様々な課題の共有と改善への協力を働きかけている。2016 年に世界で生産されている天然ゴムは年間 1,300 万トン。これは過去 20 年間で 2 倍に増加した。
今後も増加が見込まれる中、持続可能な天然ゴムの生産を拡大していくことは、イ
ンドシナ半島の環境を保全していく上で欠かせない要素となる。持続可能な天然ゴ
ムの需要と供給、二つをつなぐことで、ミャンマーの森林保全と生産者の収入向上を同時に実現させることを目指している。［3］ミャンマーは 1990 年から 2020 年にかけて、27%の森林が失われ、その中の 3 分の１は違法伐採による森林減少である。［4］2014 年、ミャンマー政府は違法伐採について、厳密に規定した。ポール・サロペックは、「政府は野生チークの輸出禁止を課した。
中国の違法伐採者 153 人を逮捕し、20 年の禁固刑を公表した。」（Paul Salopek，エ68

JSST 2022
［5］

ディター・ジャーナリスト，アメリカ）を挙げている。

森林局は、汚職防止のため、委員会の公共対応フィードバックプログラム、Public
Feedback

Program（PFP）ウェブポータルを利用して、森林局からのニュース、情報提供サービスを受けている人に、「森林局のサービスに関するパブリックコメントに対応する」ための質問を SMS でのメッセージとして送付している。また、汚職と違法伐採について、訴えがあった時、早速な対応に取り込んでいる。

4. おわりに
本稿では世界的な森林破壊の現状を示し、ミャンマーにおける森林破壊の原因として、天然ゴム栽培と違法伐採について検討した。二つの原因に関して WWWF のプロジ ェ ク ト と 環 境 保 全 林 業 省 の 公 共 対 応 フ ィ ー ド バ ッ ク プ ロ グ ラ ム 、 Public
Feedback

Program（PFP）ウェブポータルについて述べた。今後、ミャンマーの森

林破壊減少のため、違法伐採の防止、ミャンマーの森林破壊と自然環境を守る天然
ゴム栽培重視し、WWF の活動に取り込んでいく。

5. 参考文献［１］“တရ ားမ၀င်သစ်ဖမ်ားဆားရမိမ”ှု Ministry of Natural Resources andEnvironmental Conservation (Forest department) Myanmar ,http://www.forestdepartment.gov.mm （2019 年 5 月 5 日），(2022 年 7 月 12 日参
照)［2］“Myanmar Deforestation Rates & Statistics/GFW-Global Forest Watch”，
http://www.globalforestwatch.org (2022 年 6 月 13 日参照)
［3］“WWF の活動，ミャンマーで森林保全と持続可能な天然ゴム生産のプロジェク
トが始動”，WWF ジャパン http://www.wwf.or.jp/activities/activity/3839.html（2018 年 12 月 25 日），(2022 年 6 月 14 日参照)
［4］“Illegal logging and associated trade in Myanmar”, Forest Trends
Association http://www.forest-trends.org (2022 年 7 月 12 日参照)
［5］“Myanmar’s forests under pressure from illegal logging,smuggling, ” https://www.rfa.org/english/commentaries/forests-myanmar10112021104031.html/ampRFA (2022 年 7 月 12 日参照)

69

JSST 2022

Analysis of Hydrogen Bonds inside Reverse Micelles
by Molecular Dynamics Simulations
Yuri Ono1,* , Tomoko Mizuguchi2 , Mao Fukuyama3
1 Graduate School of Science and Technology, Kyoto Institution of Technology
2 Faculty of Materials Science and Engineering, Kyoto Institution of Technology
3 Institution of Multidisciplinary Research for Advanced Materials, Tohoku University∗ m2672009@edu.kit.ac.jp

The full paper of this paper has been published in Transaction of the Japan Society for
Simulation Technology.
DOI: https://doi.org/10.11308/tjsst.15.36

70

JSST 2022

コンプトンカメラを用いた放射線源位置の推定シミュレーション
Simulation of radiation source location estimation
using Compton camera
大堀 祥己*, 山口 克彦
Yoshiki Ohori*, Katsuhiko Yamaguchi

福島大学 共生システム理工学研究科
Graduate School of Symbiotic Systems Science and Technology, Fukushima University*
s2270009@ipc.fukushima-u.ac.jp

Abstract. A Compton camera is a measuring instrument that can estimate the location of radiation sources and visualize the distribution of radioactive materials. We thought that if we
could perform an analysis based on the principle of Compton camera, we could improve theaccuracy of the source distribution estimation. In this study, we attempted to estimate thesource location by analyzing radiation behavior using simulations based on particle transport
calculations.

Keywords: Compton camera, Compton cones, scattering angle

1. 諸言
除染や廃炉作業の効率化には放射線源分布の把握が必要である．線源分布を推定できる装置としてコンプトンカメラが存在する．以前までは，ピンホール式のガンマカ
メラで線源分布の推定を行っていたが，透過力の高いガンマ線を検出するために検出素子を厚くする必要があった．さらに，視野以外の方向から飛んでくるガンマ線を遮蔽する必要があり，鉛などの重い遮蔽材を使用しなければならず，重量が重くなる問題があった．しかし，コンプトンカメラはその原理から遮蔽材を必要とせず，小型化・軽量化に有利で視野が広いという特徴がある．線源位置の推定に用いられているコン
プトンカメラは分布の解析が自動で行われるものが多い．そこで，より詳細な測定や74

JSST 2022
得られたデータを応用する際には，原理に基づいた解析が有用になると考えた．本研
究では粒子輸送計算コード PHITS[1]を用いて光子の挙動解析を行い，実際のコンプト
ンカメラと同様に放射線源位置の推定が可能か検証を行った．

2. 手法
設定したシミュレーション体系について説明する．線源は 137Cs を想定し，エネル
ギー0.662 MeV の光子とした．散乱体のサイズは 1×1×3 cm3，吸収体のサイズは 1
×1×6 cm3 に設定し，それぞれ 5×5 個に並べた．散乱体を線源から 5 cm 後方に配置
し，吸収体を散乱体から 3 cm 後方に配置した．散乱体および吸収体の素材は NaI とした．設定した体系の概略図を Fig.1 に示す．計算の試行回数は 2000000 回に設定し，計算後に飛跡の様子を表す図とエネルギースペクトルのグラフが出力されるよう設定した．この時の飛跡の様子とエネルギースペクトルのグラフは Fig.2 のように示される．

Fig.1 Schematic diagram of the set system

Fig.2 Track appearance and energy spectrum in each region
本研究では，Fig.2 のように線源から放射状に出てくる 2000000 本の光子の中から，設定した条件を満たす光子に制限し，計算するようにした．例えば，Fig.1 のように
領域 112 を通り領域 213 に入射する光子に制限して計算すると，飛跡の様子とエネル
ギースペクトルは Fig.3 のように示される． 今回は領域 112 から領域 202,207,210，
211,213,214,217,222 を通る 8 通りの光子の挙動解析を行った．

75

JSST 2022

Fig.3 State of photon tracks passing thorough reg112 to reg213
and enegy spectrum in each region
散乱角は，エネルギー及び運動量保存則から算出される以下の式から求めた[2]．
𝑐𝑜𝑠𝜃 = 1 −𝑚𝑒 𝑐 2𝑚𝑒 𝑐 2(1)+𝐸2𝐸1 + 𝐸2

E1 は散乱体での検出エネルギー，E２は吸収体での検出エネルギー，mec2 は電子の静止質量エネルギーを表す．

3. 結果
解析によって得られた散乱角を Table1 に示す．
Table 1 Energy in each detector and scattering angle

今回の解析では計 8 つの散乱角を導出することができた．reg112~210 と reg112~214
では指定した領域を通る光子が存在しなかったため，散乱角を得ることはできなかった．次に，求められた散乱角から，可視化ソフト MicroAVS を用いてコンプトンコーンを作成した．

76

JSST 2022

(a)

(b)
Fig.5 Overlaid figure of 8 Compton cones(c)

得られた 8 つのコンプトンコーンを重ね合わせた図を Fig.5 に示す．Fig.5 に見られる 2 つの直線の交点は実際に線源が存在する位置を表している．コンプトンカメラの原理より，コンプトンコーンを重ね合わせた際にできる交点から線源の位置を推定す
る．8 つのコンプトンコーンを重ね合わせたところ，Fig.5(a)に示されている赤い円部
分で多く交わっていることから線源の位置は赤い円部分に存在すると推測される．実際に線源を置いた位置は 2 直線の交点である．よって，本研究では線源の存在する方向を推定できることが分かった．
今回の推定の結果，線源の正確な位置までは分からなかったが，光子の飛来してきたある程度の方向を把握することができた．推定の精度を向上させる方法として，重ね合わせるコンプトンコーンの数を増やすことが考えられる．また，推定した線源の位置と実際の線源の位置に差異が生じた原因の 1 つに，設定した散乱体および吸収体のサイズが大きいために角度分解能が低いことが挙げられる．そのため，今後の推定においては，散乱体および吸収体のサイズ設定を小さくすることを検討している．

参考文献[1] Tatsuhiko Sato, Yosuke Iwamoto, Shintaro Hashimoto, Tatsuhiko Ogawa, Takuya Furuta, Shin-ichiro Abe, Takeshi Kai, Pi-En Tsai, Norihiro Matsuda, Hiroshi Iwase, Nobuhiro Shigyo, Lembit Sihver and Koji Niita: Features of Particle and Heavy Ion
Transport code System (PHITS) version 3.02, J. Nucl. Sci. Technol, 55:6 (2018),
PP.684-690[2] Hideaki Katagiri, Wataru Satoh, Ryoji Enomoto, Ryo Wakamatsu, Takara Watanabe,
Hiroshi Muraishi, Mika kagaya, Souichi Tanaka, Kiyoto Wada, Manobu Tanaka andTomohisa Uchida: Development of an all-sky kamma-ray Compton camera based on
scintillators for high-dose environments, J. Nucl. Sci. Technol, 55:10 (2018),
PP.1172-1179

77

JSST 2022

機械学習を用いた遮蔽体分布の凝集度に対する推定精度の検証
Verification of Estimation Accuracy for Shield
Distribution Aggregation using Machine Learning近藤 勇斗 1*, 松本 正晴 2, 鈴木 健司 1, 山口 克彦 1
Yuto Kondo1*, Masaharu Matsumoto2, Kenji Suzuki1, Katsuhiko Yamaguchi1

1

1

福島大学大学院 共生システム理工学研究科
2
福島大学 情報基盤センター

Graduate School of Symbiotic Systems Science and Technology, Fukushima University
2
Information Technology Center, Fukushima University*s2270020@ipc.fukushima-u.ac.jp

Abstract. In this study, we estimated the radiation source distribution and the shield distribution from the γ-ray spectrum using machine learning, and examined how the difference in theaggregation degree of the shield distribution affects the estimation accuracy. The data formachine learning was created by Monte Carlo simulation, and the γ-ray spectrum was input
to the training data, and the distribution of the radiation source and shield was input to theanswer data. As a result, we showed that the larger the aggregation degree of the shield distribution, the higher the estimation accuracy.
Keywords: Machine Learning, Monte Carlo Simulation, Fukushima Daiichi Nuclear Power Plant
accident

1. 序論
2011 年に発生した福島第一原子力発電所の事故により，原子炉の構成物が核燃料と溶け合わさり燃料デブリを形成した．デブリは強い放射線を発するため，廃炉措置
において，人が立ち入ってデブリの位置を調べることはできない．そこで，多量な線
源分布のパターンからγ線スペクトルを計算し，学習データとしてγ線スペクトル，
解答データとして線源分布を機械学習にかけることによって，線源分布を推定する手法が提案された[1]．文献[1]では，線源と遮蔽体を設置したシミュレーションモデル
を計算し，データとして機械学習に用いている．推定対象は線源のみであり，遮蔽体の大きさを固定したデータでの推定が可能であると確認されている．本研究では，複雑なパターンの遮蔽体で推定を行うため，線源と遮蔽体の両方を推定対象とし，遮蔽

78

JSST 2022
体はランダムに分布させた．推定結果は遮蔽体分布の凝集度ごとに平均二乗誤差を取り，推定精度にどのような影響を及ぼすのか検証した．

2. 計算手法
学習データを作成するため，線源分布と遮蔽体分布を設定し，γ線スペクトルを計算した．計算には， モンテカルロ法によって 放射線の挙動をシミュ レーションする
PHITS[2]というプログラムを用いた．シミュレーションに用いたモデルは，図 1 に示
すように，1cm×1cm のメッシュで構成された 10cm×10cm の 2 次元平面を線源分布面とし，1cm 直上にγ線スペクトルの検出面を設置した．線源は 137Cs とし，メッシ
ュ中に 3 個分布させた．遮蔽体は材質を鉛とし，線源分布面と検出面との間に，10cm
×10cm の範囲に大きさ 1cm×1cm で 9 個設置した．9 個の遮蔽体は，3cm×3cm から10cm×10cm の凝集度で分布させた．図 2 は，遮蔽体分布の決定の一例である．まず，10cm 四方の分布範囲全体において，7cm 四方に遮蔽体の分布する凝集度を定める．
次にモデル中で凝集する位置をランダムに決定し，遮蔽体を 9 個分布させる．

図1

図2

シミュレーションモデルの模式図

遮蔽体分布決定の一例

線源分布と遮蔽体分布の推定を行うため，機械学習アルゴリズムの一種であるニューラルネットワークを用いた．図 3 は学習に使用したニューラルネットワークである．
1 つのパターンにつき，0.0 から 0.7 MeV までγ線スペクトルを 0.1 ごとに計 100 メッ
シュに入力するため，入力層は 700 ノード，中間層は入力層に合わせて 700 ノード，
出力層はそれぞれ 100 メッシュに収まる線源分布と遮蔽体分布と対応するため 200 ノード，計 3 層とした．
学習に用いた，遮蔽体分布範囲 4cm×4cm から 10cm×10cm のデータ，それぞれ3000 パターン計 21000 パターンを学習パターンとした．学習する際，ニューラルネ
ットワークの動作のテストに，遮蔽体分布範囲 10cm×10cm のデータ 3000 パターンを用い，テストパターンとした．
線源分布と遮蔽体分布の推定のため， 3cm×3cm から 10cm×10cm 遮蔽体分布範
囲のパターンをそれぞれ 5 パターン入力した．これを推定パターンとした．
なお，テストパターンと推定パターンは，学習パターンには含まれていない．

79

JSST 2022

図 3 学習に用いたニューラルネットワーク

3. 結果
ニューラルネットワークの学習回数に対する誤差を図 4 に示した．学習パターンにおける誤差を学習誤差，テストパターンにおける誤差をテスト誤差とした．図 5 は，遮蔽体分布の凝集する範囲 3cm×3cm で計算したγ線スペクトルから推定した線源分布と遮蔽体分布である．図 6 には，遮蔽体の個数を分布範囲で割った凝集度に対して，遮蔽体分布の解答と推定結果の平均二乗誤差を取ったグラフを示した．

図4

ニューラルネットワークの学習誤差80

図 5 線源分布と遮蔽体分布の推定結果

JSST 2022

図6

遮蔽体分布の平均二乗誤差

図 4 より，学習データの誤差に比べて，テストパターンにおける誤差に変動が見られなかった．これはテストパターンの数が学習パターンに比べて非常に少なかったためと考えられる．図 5 からは，遮蔽体分布の推定精度が線源分布よりも低いことが目視で確認できる．遮蔽体の位置は，γ線が遮蔽体に衝突し散乱することによって推定されるため，線源が少ないためγ線の本数も少なく，解答にはない位置に遮蔽体を推定した可能性がある．
遮蔽体分布の平均二乗誤差は，分布の凝集度が大きくなるにしたがい減少していることが図 6 から読み取れる．したがって，凝集度が大きいほど遮蔽体の推定精度は高いと言える．
今後の課題として，推定パターンは遮蔽体分布範囲ごとに 5 つと少ないため，100
パターンぐらいまで増やして凝集度に対する誤差を算出したい．加えて，鉛と異なる材質の遮蔽体でモンテカルロ計算を行ったデータで分布を推定し，凝集度に対する推定精度の変化を検証したい．

参考文献[1] Takero Uemura, Katsuhiko Yamaguchi: Estimation of radiation source distribution using machine learning with γ ray en-ergy spectra, Journal of Advanced Simulation in
Science and Engineering, 7:1 (2019), 71-81.
[2] Tatsuhiko Sato, Yosuke Iwamoto, Shintaro Hashimoto, Tatsuhiko Ogawa,Takuya Furuta, Shin-ichiro Abe, Takeshi Kai, Pi-En Tsai, Norihiro Matsuda,
Hiroshi Iwase, Nobuhiro Shigyo, Lembit Sihver, Koji Niita: Features of Particle
and Heavy Ion Transport code System (PHITS) version 3.02, Journal of Nuclear
Science and Technology, 55:6 (2018), 684-690.

81

JSST 2022

深層学習を用いた静電場シミュレーション

Electrostatic ﬁeld simulation using deep learning
相良 翼 1,* , 松本 正晴 2 , 鈴木 健司 1 , 山口 克彦 1
Tsubasa Sagara1,* , Masaharu Matsumoto2 , Kenji Suzuki1 , Katsuhiko Yamaguchi1
1 福島大学 共生システム理工学類
2 福島大学 情報基盤センター
1 Faculty of Symbiotic Systems Science, Fukushima University
2 Information Technology Center, Fukushima University∗ s1910077@ipc.fukushima-u.ac.jp

本論文のフルペーパーは，日本シミュレーション学会論文誌に掲載されています．

DOI: https://doi.org/10.11308/tjsst.15.20

82

JSST 2022

メゾスコピック領域における磁性体のマルチスケール解析
Multiscale analysis of magnetic materials in the mesoscopic region
五十嵐靖幸*, 山口克彦
Yasuyuki Igarashi*, Katsuhiko Yamaguchi
福島大学大学院 共生システム理工学研究科
Graduate School of Symbiotic Systems Science and Technology, Fukushima University*

s2170005@ipc.fukushima-u.ac.jp

Abstract. In magnetization analysis, the crystal structure can be considered by the Monte
Carlo method. On the other hand, analysis using the LLG equation makes it easy to calculatemagnetic clusters of a size that is difficult to calculate with the Monte Carlo method. In orderto analyze the magnetization of the magnetic material in the mesoscopic region considering
the crystal structure by combining the two analysis methods, the magnetization processes ofeach method were compared.

Keywords: micro magnetic simulation, LLG equation, Monte Carlo method

1. 諸言
磁化過程解析をおこなう方法として磁化の時間変化を計算する際に用いられる
Landau-Lifshitz-Gilbert 方程式(LLG 方程式)を使用するものがある. LLG 方程式を用いた解析では磁性体を大きさが一定かつ磁化が一様なセルで分割し, セル 1 つに対して1 スピンとして計算を行う. また, 他の解析方法としてはモンテカルロ法が挙げられる. モンテカルロ法では原子 1 つを 1 スピンとして, 乱数を用いて計算を行う. 原子 1
つを 1 スピンとして計算を行うことから結晶構造を考慮して計算が行える. しかし,磁性体の大きさが大きくなると計算するスピン数が多くなり, 時間が膨大にかかり,容易に計算が行えなくなってしまう. これよりモンテカルロ法はミクロ領域の計算が有効といえる. 一方, LLG 方程式を用いた解析では, 結晶構造を考慮せず 1 つのセルを 1 スピンとみなして計算を行うので, 同じ大きさの磁性体を想定した場合, モン
テカルロ法よりも計算が行いやすくなる. また, より大きなクラスタについても計算が可能なことから, LLG 方程式を用いた解析ではミクロ領域より大きなメゾスコピッ
ク領域での計算も有効といえる.
よって, モンテカルロ法と LLG 方程式を用いることで, 結晶構造を考慮しつつメ

86

JSST 2022
ゾスコピック領域における計算が可能なマルチスケールでの磁化解析が行えるのではないかと考えられる. 本研究では, モンテカルロ法と LLG 方程式を用いて, それぞれについて解析を行い磁化過程の結果を比較する.

2. 手法
まず LLG 方程式について(1)式に示す.
𝑑𝑴⁄𝑑𝑡 = −𝛾𝑴 × 𝑯𝒆𝒇𝒇 + (𝛼 ⁄𝑀𝑠 )𝑴 × 𝑑𝑴⁄𝑑𝑡(1)(1)式中の𝑴は磁化ベクトル, γはジャイロ磁気定数(𝛾 = 1.76 × 107 [1/𝑂𝑒 ⋅ 𝑠]), 𝑯𝒆𝒇𝒇 は有効磁場, αはギルバート減衰定数(0 ≤ 𝛼 ≤ 1), 𝑀𝑠 は自発磁化である. 𝑯𝒆𝒇𝒇 は,磁性体中で隣接する磁気モーメントに働く交換相互作用による交換磁場, 結晶磁気異方性による異方性磁場, 磁性体中の磁気モーメント同士に働く磁気双極子相互作用による磁気双極子磁場, 外部磁場の影響を考慮した.
モンテカルロ法は乱数を用いて何度も試行を繰り返し, 近似解を求める方法である. 本研究ではその中でも, メトロポリス法を用いて解析を行った. メトロポリス法
は, 想定する磁性体クラスタ内のスピンの向きをランダムに変えて, 向きを変えた前
後のハミルトニアンが小さい場合のスピンの向きを採用し, それを繰り返す手法である[1]. ハミルトニアンは LLG 方程式の有効磁場と対応させて, 交換相互作用による交換エネルギー, 結晶磁気異方性による異方性エネルギー, 磁気双極子相互作用による磁気双極子エネルギー, 外部磁場によるゼーマンエネルギーを考慮した.

3. 結果
解析には LLG 方程式を用いた計算, モンテカルロ法での計算それぞれについて C
言語で作成したプログラムを使用した. LLG 方程式を用いた解析では, サイズが40×2×2[nm3]の直方体クラスタを想定して, 分割するセルの大きさは 2×2×2[nm3]で行った. また, モンテカルロ法による解析では BCC 構造を考慮して 10×0.5×0.5[nm3]の直方体クラスタを想定した. それぞれの解析において材料は鉄を想定した. 外部磁場
は-1,000,000[Oe]～1,000,000[Oe]の間で変化させて印加し磁化反転させて磁化過程を調べた. LLG 方程式を用いた計算による磁化曲線を Fig.1, モンテカルロ法を用いた計算による磁化曲線を Fig.2 に示す.

Fig.1 Magnetization curve by LLG equation

Fig.2 Magnetization curve by Metoropolis method

87

JSST 2022
図より, LLG 方程式を用いた計算とモンテカルロ法を用いた計算それぞれにおいて磁化反転を伴う磁化過程が確認できた. 2 つの方法において, それぞれ磁化が安定していて計算が比較的容易な磁性体クラスタで解析を行った. 今後パラメータや想定するクラスタの形状を調整することで, 2 つの解析方法において同じクラスタの条件
で磁化過程の結果が一致することが期待でできる. また, 結晶構造の乱れと考慮した計算を行い, 磁化過程だけでなくエネルギーについても評価を行うことで, 鉄鋼材料の劣化具合についての評価に繋げられることが期待できる.

参考文献[1] 三俣千春, 磁気工学の解析法, 共立出版, pp.119-121, 2013.

88

JSST 2022

サロゲートモデルによるデータ拡張の効用について

The Benefit of Data Augmentation Performed by Surrogate Modeling野村拓海
Takumi Nomura

11

東京大学大学院 工学系研究科 システム創成学専攻

Department of Systems Innovation, School of Engineering, The University of Tokyo
nomura-takumi11026@g.ecc.u-toyko.ac.jp

Abstract. This study proposes a method of data expansion in the analysis of inverse problems.
An inverse problem is a problem in which the cause is estimated from the measurement datadetected, as opposed to a forward problem in which the result is estimated from the cause orcondition. There are many examples of inverse problems in the field of materials and structures,and machine learning has been used for their analysis in recent years. For example, machinelearning has been used in various cases such as nondestructive inspection of structures[1] andidentification of impact loading conditions[2].
Since inverse problems are often estimated with insufficient parameters, they require a larger
number of data for analysis than simple forward problems. For this reason, physical simulations
are often used to prepare the data set.
However, the computational cost of physical simulation is high, and this can be a major obstacle
when large-scale analysis is required within a limited time frame or when high-performancecomputers are not available.
As a solution to this problem, this study focuses on the method proposed by G. R. Liu (2019),
which uses surrogate model data expansion as an alternative to physical simulation [3]. In thisstudy, we will test the effectiveness of this proposed method using data from structural internal
defect estimation.
Translated with www.DeepL.com/Translator (free version)
Keywords: Surrogate Model, Machine Learning, Inverse Problem, non-Destructive Inspection

1. 序論
本研究は逆問題の解析におけるデータ拡張の手法を提案するものである。

89

JSST 2022
逆問題とは原因や条件からその結果を推定する順問題に対し、逆に検出される測定データからその原因を推定する問題のことをいう。材料・構造物分野では逆問題を扱う事例が数多く存在し、近年その解析に機械学習が用いられることが多い。例えば構造物の非破壊検査[1]や衝撃荷重条件の同定[2]など様々な事例において成果を出している。
逆問題はパラメータが不足した状態で推定を行うことが多いため、単純な順問題の解析に比べて解析に多くのデータ数を必要とする。そのため物理シミュレーションによりデータセットを揃える方法が取られることが多い。
しかし、物理シミュレーションによる計算コストは大きく、大規模かつ限られた時間内での解析が必要な場合や、高性能な計算機が揃えられない場合に、これが大きな障壁となり得る。
これに対する解決策として本研究では、物理シミュレーションの代替としてサロゲートモデルによるデータ拡張を行うという G. R. Liu(2019)が提案している手法に着目する[3]。今回は構造物の内部欠陥推定のデータを用いてこの提案手法の効果を検証する。
キーワード：サロゲートモデル、機械学習、逆問題、非破壊検査

2. 既存研究の手法紹介
2.1

サロゲートモデルとはサロゲートモデルとは、代替モデルとして物理現象の解析を模擬するモデルのことをいう。今回は物理現象をシミュレーションで解析する代わりに、ニューラルネットワーク
を使ってこの過程を代替するモデルを作りデータ生成に用いる。ここでサロゲートモデ
ルは、原因からその結果が生じる物理現象を模しているため順問題を解く順解析モデル
である。

2.2

順解析モデルの逆解析モデルへの適用

順解析モデルを逆解析モデルの学習に利用する方法としては 2 つのアプローチが考えられる。
一つは順解析モデルの重みやバイアスを利用して逆解析モデルで直接計算を行うアプ
ローチである。これは実質ネットワークが繋がって一つになるため計算コストの削減とリアルタイムでの順解析・逆解析の並列化が期待できるアプローチである。しかし、ニューラルネットワークの隠れ層の活性化関数には通常非線形関数を使用するものの、重みやバイアスの反映は非線形性に弱く、どうしても逆解析モデルの再定式化が必要になってしまうため実装は難しい。
もう一つのアプローチは、サンプルを作成するサロゲートモデルとして利用する順解析ネットワークと、そのサンプルを訓練データとして学習させる逆解析ネットワークとを別々に組むものである。順解析と逆解析とで別々にネットワークを実装することは難しくない。本研究ではこちらのアプローチを採用する。

3. 機械学習モデル構築
3.1

サロゲートモデルの構築

サロゲートモデルの学習は、入力を欠陥の大きさと位置、出力を構造物表面の温度分布

90

JSST 2022とした回帰問題として行った。
モデルの構造としては入力層を 2 つ用意して結合し、全結合層を多層に重ねて出力層に繋がるモデルとした。

3.2

逆解析モデルの構築

逆解析モデルの学習は、入力を構造物表面の温度分布データ、出力を欠陥の大きさ、欠陥の位置とし、大きさ推定と位置推定は別々の学習行った。いずれも分類問題として扱った。
モデルの構造としては、分類問題に適した畳み込みニューラルネットワークを用いた。
入力層から畳み込み層、プーリング層での圧縮を経て、全結合層、出力層へと繋がるモデ
ルとした。

4. 実験
4.1

データの概要今回使用するデータは構造物内部の欠陥を模した 1m 四方の平面の全位置について大きさ 0.02〜
0.1m 四方の欠陥を 1 つ設定したデータと、これに有限要素法による温度分布シミュレーションを施
図 1 温度分布データ例
すことで作成した図 1 のような温度分布データである。データの総数は 10135 個であり、本研究ではそのうちの一部を抽出し小さいデー
タ数で学習を行う。
具体的には、
1000,2000,3000,4000,5000 個の抽出をそれぞれ行った。
値として用いるのは欠陥の大きさ、欠陥の位置、構造物の表面温度分布の 3 つのデータである。
欠陥の大きさは 0.02,0.04.0.06.0.08,0.1m のデータを使用する。欠陥の位置はサロゲー
トモデルでは座標データを、逆解析モデルではある程度の幅を持たせた推定の正答率を測るため位置を 5×5=25 個の範囲に分類したデータを用いた。表面温度分布は図 1 のように構造物表面を表す左辺 1 列 50 個の温度データを用い、その際より特徴が現れるよ
う、欠陥のない構造物の温度分布データとの差を学習に使う値として採用した。

4.2

実験 1：使用データ数によるデータ拡張の効果検証

表 1 逆解析モデル評価比較(データ数による比較)

実験 1 では、逆解析モデルを用いて、使用データ数によるデータ拡張の効果の違いを、大きさ推定と位置推定のそれぞれについて検証した。拡張倍率は 2 倍で統一し、データ抽出数は前述に加え 250 個、500 個の 2 つの場合を追加した。

91

JSST 2022
その結果、大きさ推定ではいずれのデータ数の場合も十分に精度が向上した。一方位置推定では使用データ数 2000 個以上の時はむしろ精度が下がり、1000 個以下の時に精度が十分に向上した。

4.3

実験 2：データ拡張倍率によるデータ拡張の効果検証

表 2 逆解析モデル評価(データ拡張倍率による比較)
実験 2 では、データ拡張倍率によるデータ拡張の効果の違いを検証するため、2 倍拡張に加えてデータを最大の 10135 個に拡張した場合(ここでは max 拡張時とする)との比較を行った。
その結果、大きさ推定ではいずれのデータ数の場合も max 拡張時の方が 2 倍拡張時よりも精度が高くなった、一方、位置推定では使用データ数が 500 個と 1000 個の場合に 2倍拡張時の精度が max 拡張時の精度を上回った。拡張倍率が大きいほどデータ拡張の効
果が増すわけでは必ずしもないということがわかる。
また、実験 1、2 のいずれについても大きさ推定と位置推定で結果に大きな差が生まれているが、これは各問題の性質やハイパーパラメータの設定の差に起因していると考え
られる。

5. 結論
上記実験から以下のような結論を得た。
・小さいデータ数でのモデル評価は可能であるが値を外す場面も見受けられた・データ拡張はある程度小さいデータ数において効果を発揮し、同数の拡張前デー
タを使用した場合と比べても遜色のない精度を出すことが できた
・拡張倍率を大幅に高めることは精度を上げることも下げることもあり、必ずしも高い拡張倍率が好ましいわけではない
・入力と出力の関係や問題の性質によってデータ拡張の効果に差が生まれる
今後の課題としては、シミュレーションデータの前処理、回帰問題としての評価、
ハイパーパラメータ最適化等によるモデル選択方法の検討、が挙げられる。これらの解決によりデータ拡張の効用についての更なる一般化を進めていく。

参考文献[1] 阪上隆英, 溶接学会誌, 72(4), pp. 251-255, 2003.
[2] 中村正行, 福島斉春, 日本機械学会北陸信越支部総会・講演会講演論文集, 45, pp.
447-448, 2008.
[3] G. R. Liu, International Journal of Computational Methods, 16(8), 2019.

92

JSST 2022

単成分ガラス形成系における配置エントロピー計算
Configurational Entropy of a Monatomic Glassy System
上野 綾大 1*, 水口 朋子 2Ayata Ueno1*, Tomoko Mizuguchi2

1

京都工芸繊維大学 大学院工芸科学研究科 材料制御化学専攻
2
京都工芸繊維大学 材料化学系

Graduate School of Science and Technology, Kyoto Inssitute of Technology
Faculty of Materials Science and Engineering, Kyoto Institute of Technology* m1672004@edu.kit.ac.jp

12

Abstract. The glass transition is one of the major unsolved problems in condensed matterphysics. A theoretical framework that can explain the glass transition in a unified manner
does not yet exist, and various theories are currently being proposed. One of the most promising candidates is the theory based on the configurational entropy, which is the entropy corresponding to the number of atomic configurations. Various methods have been proposed toestimate the configurational entropy in computer simulations. In general, it is not easy to calculate only the configurational entropy, since the total entropy is composed of multiple contributions including configurational, mixing, harmonic vibrational, and non-harmonic vibrational components. Therefore, the relationship between the configurational entropy and thestructure is still a central issue for the study of glass transition. To discuss this relationshipquantitatively, we calculate the configurational entropy in simulations using a monatomic
glass-forming model where the atoms interact each other isotropically.
Keywords: Molecualr dynamics simulations, Configurational entropy, glass, Adiabatic switching,Reversible scaling, Isotropic monatomic system

1. 序論
ガラス転移とは、液体を急冷することで結晶化をしそこない、ランダムな分子配置で凍結してしまう現象である。ガラス転移を統一的に説明できる理論的枠組みはまだ存在しておらず、現在様々な理論が提唱されている。そのうち、原子配置の数に対応する配置エントロピーを用いた理論は有力候補の 1 つであり、様々な配置エントロピーの計算手法が提案されてきた。一般に、全エントロピーは、配置エントロピー、混合エントロ
ピー、調和振動エントロピー、非調和振動エントロピーなど、複数の寄与から構成されるため、配置エントロピーだけを計算することは容易ではない。ガラス形成モデルは通常、ガラス構造の安定性の観点から多成分系が用いられる。今まで 2 成分 Lennard-Jones
ポテンシャル（Kob-Andersen モデル）や合金用の EAM（Embedded-Atom Method）

93

JSST 2022
ポテンシャルにおける配置エントロピー計算の結果[1,2]が報告されているが、いずれのモデルも多成分系であり混合エントロピーの寄与を取り除く必要がある。本研究では、
配置エントロピーに対する知見の蓄積のため、これまで計算されたことが無い単成分モデルを用いて配置エントロピーを計算し、ガラス転移と構造の定量的な関係を議論する。

2. 計算手法
等方的に相互作用する単原子ガラス形成系である Lennard-Jones-Gauss（LJG）ポテンシャル[3]を用いて分子動力学（MD）シミュレーションを行い、配置エントロピーを計算する。ガラス形成物質である合金の有効 2 体ポテンシャルの一般形は、近距離
斥力項と減衰振動項（Friedel 振動項）からなる振動型ポテンシャルである。LJG ポテンシャルはこのよ
うな振動型ポテンシャルを第 2 極大でカットオフしたものとみなすことができる（Fig. 1）。以下、物理Fig. 1. LJG ポテンシャルの概形.
量は全て Lennard-Jones 単位系で無次元化された値
横軸は２粒子間の距離、縦軸はポを用いる。
テンシャルエネルギーを表し、ど広い温度領域での自由エネルギーを計算するこ
ちらも無次元化されている.
とで配置エントロピーを求める。単原子モデルにおいては混合エントロピーの寄与を無視できるため、
配置エントロピーは全エントロピーから振動エントロピー（調和振動エントロピーおよび非調和振動エントロピー）を引いたものだと定められる。配置エントロピーの計算には、Adiabatic Switching（AS）法[4]と Reversible Scaling（RS）法[5]を用いる。これらは不可逆過程における計算から配置エントロピーを求める手法である。温度 T = 0.1 でのガラスを参照系とした RS 法によって、任意の温度におけるガラスの自由エネルギーを得ることができ、この自由エネルギーの温度に関する微分を振動エントロピーとみなす。一方、
温度 T = 2.5 での液体を参照系とした RS 法によって、任意の温度における液体の自由エネルギーを得ることができ、この自由エネルギーの温度に関する微分を全エントロピーとみなす。RS 法において基準となる T = 0.1（ガラス）および T = 2.5（液体）での自由エ
ネルギーは、それぞれ Einstein crystal および Uhlenbeck-Ford model を参照系とした AS 法によって得ることができる。MD 計算には汎用古典 MD アプリケーションである
LAMMPS（Large-scale Atomic/Molecular Massively Parallel Simulator）を使用する。

3. 結果
T = 2.5 の液体から T = 0.1 まで急冷し、構造緩和させたあとの静的構造因子 S(q)を Fig. 2に示す。S(q)は原子構造の規則性を調べる上で重要な物理量である。結晶であれば S(q) は離散的なピークを有し、２次元の回折像にはスポットが現れるが、Fig.2 にスポットは見られず、ランダム構造に特有のリング形状が現れている。このことから T = 0.1 における構造は結晶ではなく、液体のようなランダム配置を持っていることがわかる。さらに、冷却中のポテンシャルエネルギーの変化を Fig. 3 に示す。エネルギーに不連続な跳びはな
く、液体から結晶への一次相転移は起こっていないことが分かる。T = 1.5 を超えたあた
りで二次相転移的な傾きの変化が見られることから、S(q)と合わせて判断すると、低温で94

JSST 2022

Fig. 3. 圧力 P = 10 におけるポテン
シャルエネルギーの温度変化

Fig. 2. T = 0.1 における構造因子 S(q)ガラス状態が形成されており、T = 1.5 付近で過冷却液体へと転移したと考えられる。T =0.1 の配置をガラスの参照構造、T = 2.5 の配置を液体の参照構造とし、AS 法および RS
法を用いて配置エントロピーを計算した結果を、講演では示す。

参考文献
[1]
[2]

[3]
[4]
[5]

L. Berthier, M. Ozawa, and C. Scalliet: Configurational Entropy of glass-forming liquids, J.
Chem. Phys. 150 (2019), 160902.
R. Alvarez-Donado and A. Antonelli: Splitting up entropy into vibrational and configurational contributions in bulk metallic glasses: A thermodynamic approach, Phys. Rev. Res. 2(2020), 013202.
T. Mizuguchi and T. Odagaki: Glass formation and crystallization of a simple monatomic
liquid, Phys. Rev. E 79 (2009), 051501.
M. de Koning, A. Antonelli, and S. Yip: Optimized Free-Energy Evaluation Using a Single
Reversible-Scaling Simulation, Phys. Rev. Lett. 83:20 (1999), 3973.
M. de Koning and A. Antonelli: Adiabatic switching applied to realistic crystalline solids:
Vacancy-formation free energy in copper, Phys. Rev. B 55:2 (1997), 735.

95

JSST 2022

Numerical analysis of hydrodynamic equations using
Spreadsheet software
Sakuya Inoue 1*, Kiyomi Yamane 2, Munkhbileg 2, Hiroki Homma 2

1

Department of Production and Construction Systems Engineering, Matsue National College of2

Technology
Department of Mechanical Engineering, Matsue National College of Technology
P2105@matsue-ct.jp

Abstract. In this study hydrodynamic equations, the continuity equation and the Navier-Stokes
equation were analyzed using a spreadsheet functions of Microsoft Excel, with an iteration by
the Newton-Raphson method. Rayleigh flow problem was solved by proposed method as an
example. And it matched well with the theoretical solution.
Keywords: Numerical analysis, Microsoft Excel, Spreadsheet function , Navier Stokes equation,
Newton’s method

1. Introduction
The purpose of this research is to establish a method for using spreadsheet software to analyze differential equations[1]. This method does not require any programming, and differential
equations can be analyzed by simply entering mathematical expressions into the cells of Excel. The method is also capable of implicit analysis, which allows a wider range of applicable
equations and shorter computation time than the explicit method. In this report, the proposed
method is applied to the Navier-Stokes equations of fluid dynamics to analyze the Rayleighproblem, and its validity is verified.

2. Numerical analysis of differential equations by Microsoft Excel
Numerical analysis of differential equations using Excel can be roughly classified into two
types: one is to write a program using VBA (Visual Basic Applications) and the other is to use a96

JSST 2022
spreadsheet function.The impilict method[1]

has come up as a new method in place of the ex-pilict method[2]. This method discretizes differential equations, applies the Newton-Raphson
method to obtain an iterative formula to solve the obtained difference equation by the implicit
method, inputs the formula into Excel cell, and finds the solution by iterative calculations.
This method can be analyzed in the same procedure regardless of the type of differential
equation. It has been applied to ordinary differential equations and partial differential equations. This study aims to apply the method to a more complex partial differential equationsuch as Navier-Stokes equation of fluid dynamics, and presents analytical results for the Rayleigh problem.

3. Unsteady analysis of Rayleigh problems
The Rayleigh problem[3] is an unsteady flow generated in the fluid near the bottom surface
when the tank is filled with fluid and the bottom surface suddenly moves from rest position tox direction at a constant velocity U, as shown in Figure 1. As shown in the figure, the flow
velocity occurs parallel to the x-axis and its value decreases as it moves away from the wall
surface which is y=0.

y

u(y)
0

0
Fig. 1Ux

Rayleigh flow problem.

This problem is solved using the conservative form of the following Navier-Stokes equation(x component), which is a fundamental equation in fluid mechanics.
𝜕𝑢 𝜕𝜕
1 𝜕𝑝 𝜕𝜕𝑢𝜕𝜕𝑢(𝑢 ∙ 𝑢) +(𝑣 ∙ 𝑢) ++−(𝜈 ) −(𝜈 ) = 0𝜕𝑡 𝜕𝑥
𝜕𝑦𝜌 𝜕𝑥 𝜕𝑥 𝜕𝑥
𝜕𝑦 𝜕𝑦

(1)

The flow velocity is u and v in the x and y directions, time is t, pressure is p, density is ρ, and kinematic viscosity is ν (=μ/ρ). For the present analysis, it is assumed that the flow does not occur
in the y direction (v=0) and that no pressure change occurs (∂p/∂x=0), and only the differential equation (1) is analyzed by the following procedure.
1) Conversion to differential equation
2) Application of the Newton-Raphson method
3) Input data to Excel and perform analysis

97

JSST 2022

3.1 Conversion to difference equation
In numerical calculations, the domain is divided into a grid. For each grid point (i, j), Eq. (1)is integrated over the neighboring rectangular domain and discretized by the finite volume
method to yield the following equation.
𝜌∆𝑥∆𝑦 𝑛+1
𝜇∆𝑦 𝑛+1
𝑛
𝑛+1 2
𝑛+1 2
𝑛+1
𝑛+1
𝑛+1
− 𝑢𝑖−1,𝑗)−− 𝑢𝑖−1,𝑗− 𝑢𝑖,𝑗+ 𝑢𝑖−2,𝑗(𝑢𝑖,𝑗 − 𝑢𝑖,𝑗) + 𝜌∆𝑦 (𝑢𝑖,𝑗(𝑢)∆𝑡
2∆𝑥 𝑖+1,𝑗
𝑛+1
𝑛+1
𝑛+1
+ 𝑢𝑖,𝑗+1− 2 ∗ 𝑢𝑖−1,𝑗𝜇∆𝑥 𝑢𝑖−1,𝑗+1
−
(
𝑛+1
𝑛+1
𝑛+1 ) = 𝐹𝑖,𝑗2∆𝑦 −2 ∗ 𝑢𝑖,𝑗+ 𝑢𝑖−1,𝑗−1+ 𝑢𝑖,𝑗−1(2)

Equation (2) becomes 0 once the solution ui, j is determined, but in mid process step, the left
side becomes Fi,j.

3.2 Application of the Newton-Raphson method
This method finds the solution by decreasing the residual Fi, j through iterative calculations.When applied to equation (2), the new velocity ui, jn+1 at time n+1 can be obtained. The
principle equation is as follows, assuming the number of iterations to be k.

𝑘+1 𝑛+1
𝑛+1
𝑢𝑖,𝑗 = 𝑘+1𝑢𝑖,𝑗
−

𝐹𝑖,𝑗
𝜕𝐹𝑖,𝑗
𝑛+1
𝜕𝑢𝑖,𝑙

(3)

Substituting the left-hand side of equation (2) for the residuals Fi, j in the principle equation (3)yields the following equation.

𝜌∆𝑥∆𝑦 𝑛+1𝑛∆𝑡 (𝑢𝑖,𝑗 − 𝑢𝑖,𝑗 )
𝑛+1 2
𝑛+1 2
+𝜌∆𝑦 (𝑢𝑖,𝑗− 𝑢𝑖−1,𝑗)𝜇∆𝑦 𝑛+1
𝑛+1
𝑛+1
𝑛+1− 𝑢𝑖−1,𝑗− 𝑢𝑖,𝑗+ 𝑢𝑖−2,𝑗)2∆𝑥 (𝑢𝑖+1,𝑗
𝑛+1
𝑛+1
𝑛+1
𝜇∆𝑥 𝑢𝑖−1,𝑗+1 + 𝑢𝑖,𝑗+1 − 2 ∗ 𝑢𝑖−1,𝑗− 2∆𝑦 (
𝑛+1
𝑛+1
𝑛+1 )
−2 ∗ 𝑢𝑖,𝑗+ 𝑢𝑖−1,𝑗−1+ 𝑢𝑖,𝑗−1
𝑘+1 𝑛+1
𝑘 𝑛+1
𝑢𝑖,𝑗 = 𝑢𝑖,𝑗 −𝜌∆𝑥∆𝑦𝜇∆𝑦 𝜇∆𝑥𝑛+1∆𝑡 + 2𝜌∆𝑦𝑢𝑖,𝑗 + 2∆𝑥 + ∆𝑦−

(4)

Converge equation (4) with Excel iterations to obtain the solution ui, jn+1 for all grid points (i, j).

98

JSST 2022

3.3 Enter data into Excel and perform analysis
The analysis is performed by entering equation (4) into all cells in Excel. The analyticalconditions are as follows: height h=4×10-4 [m], number of divisions Ny=40, ∆y=10-5 [m],∆t=10-3 [s] for space and time increments, the fluid used is water (density ρ=103 [kg/m3],
viscosity μ=10-3 [Pa-s]), and velocity U=1 [m/s] for the bottom surface. The boundary conditions for the lower surfaces were the adhesion conditions (u=U at y=0).
Figure 2 compares the analytical result u at t=0.01 (n=10) with the theoretical solution. In
the figure, the symbols marked with ● are the results obtained by Excel, and the solid line is
the theoretical solution (Rayleigh solution (4)). The figures alings with each other.

0.0004y [m]: Excel: Theoretical solution

0.000200
Fig. 2u [m/s]1

Unsteady flow velocity profiles predicted by the Navier-Stokes equation (1) andanalytical solution.

4. Conclusion
The method of using spreadsheet software to analyze differential equations was applied to
the Navier-Stokes equations to analyze the Rayleigh problem, and the validity of the results
was confirmed.

References[1] Yuasa, K. and Yamane, K., Numerical analysis of differential equations utilizing Microsoft Excel (Analysis method using a spreadsheet function), Proceedings of the JSME
Chugoku-Shikoku Branch Conference (2017), 1119 (in Japanese).
[2] Iwai, H., Ohmura, T., Koito, Y., Kobayashi, K., Tomimura, T., Haneda, M., Hirasawa,
S. and Yoshida, H., Thermal fluid simulation with Excel and a mouse (3nd Edition)(2010), Maruzen (in Japanese).
[3] Toda, M. "30 Lectures on Fluid Mechanics", p. 178, Asakura Shoten, (2001) (in Japanese).

99

JSST 2022

エンジン筒内流動におけるサイクル変動成分の解析
Analysis of cycle-to-cycle variation components
in engine cylinder flow
森川 静 1, 伊藤裕一 2*
Shizuka MORIKAWA1, Yuichi ITOH2*

1

木更津工業高等専門学校 専攻科 機械・電子システム工学専攻
2
木更津工業高等専門学校 機械工学科
1
Advanced Course of Mechanical and Electrical Engineering, National Institude of Technology,
2

Kisarazu College
Department of Mechanical engineering, National Institude of Technology, Kisarazu College*itoh@m.kisarazu.ac.jp

Abstract. In recent years, cycle-to-cycle variation has attracted attention in the developmentof internal combustion engines for automobiles. The cycle-to-cycle variation occurs when the
combustion state changes. This cycle-to-cycle variation can be reduced the exhaust gas composition and higher thermal efficiency. If this fluctuation can be suppressed, harmful exhaustgas can be reduced and high thermal efficiency can be obtained. Factors that cause cycle-tocycle variation include in-cylinder flow distribution, fuel concentration distribution, and incylinder temperature distribution. Therefore, in this study, we focused on the cycle fluctuationof the engine cylinder flow and analyzed the difference of the engine cylinder flow dependingon the length of the intake pipe using the 3D fluid simulation software HINOCA. As a result,
we were able to extract and evaluate the cycle-to-cycle variation component by using 3D fluid
simulation.
Keywords: Numerical analysis, internal combustion engine, cycle-to-cycle variation, engine incylinder flow

1. まえがき
現在，地球温暖化を主とする環境問題が注目されている．我々の身近にある自動車
もこの問題に大きく影響を受けており，環境に関する規制がより一層厳しくなっている．そのため，クリーンな自動車開発が求められている．エンジン開発は年々高度化しており，従来は平均流を評価していたが，近年では乱流も評価して開発を行っている．現在では平均流や乱流に加え，サイクル変動も注目されている．エンジンはその特性上，サイクル毎に流れ場や燃焼状態が変化し，出力や排ガス成分もサイクル毎に100

JSST 2022
変化する[1]．これをサイクル変動という．サイクル変動を小さくすることができれば燃焼状態の安定化による排ガスの改善や出力のばらつきが少なくなるため，環境と性能ともに改善することができる．そのため，近年ではサイクル変動が注目され，それによる影響などが研究されている．本研究ではエンジンの筒内流動に注目し，数値解
析によるシミュレーションを用いてサイクル変動を解析する．

2. 解析方法
2.1. 使用ソフトウェア
本研究で用いるシミュレーションソフトウェアとして HINOCA[2]を用いる．
HINOCA の特徴は，直交格子法及び境界埋め込み法に基づいた，非定常圧縮性 3 次元
流体解析ソルバであり，メッシュ生成を必要としないため，短時間で計算を投入することができることである．燃焼版と流動版があり，本研究では流動版を用いた．

2.2. 解析モデル
エンジンモデルは 4 バルブ 4 ストロークエンジンであり，シリンダ，ピストン，吸排気バルブ，吸排気管を組み合わせたものである．図 1 に解析モデルを示す．吸気時の乱れた流れがエンジン筒内の流れに与える影響を調べるため，吸気管を延長し，乱流促進体を設けた．吸気管の有無による比較として，延長吸気管を付けていないモデル(0D)と，吸気管直径の 4 倍の長さの延長吸気管を取り付けたモデル(4D)との比較を行った．図 2 に解析モデルの外観を示す．
Valve

ExhaustInflow
CylinderFig.2. Comparison of analysis models(4D,0D)

Fig.1. Analysis model

2.3. 解析条件
表 1 に解析条件を示す．また，解析モデル内のモニタ点での時系列データを取得した．モニタ点はエンジン筒内の点火プラグ直下の位置に設定した．

101

JSST 2022
Condition name
Cell size [mm]Number of cells
Turbulence modelTime integration method
CFL
Angle step [deg]IntakeBoundary condition

Wall
ExhaustTable 1. Analysis condition
0D

4D
1.0

3,328,000

3,428,352
WALE
LU-SGSca. 30ca. 3.0
0.5
0.05
Total pressure and tempera- Total pressure and temperature change inflow
ture fixed inflow
Isothermal non-slip
Static pressure fixed outflow

3. 解析結果
流速の可視化結果を図 3 に示す．図は 4 サイクル目の吸気行程時のバルブ中央断面を示している．これより，エンジン筒内の流動分布が 4D モデルと 0D モデルで変化していることがわかる．これは乱流促進体によって発生した乱流分布がエンジン筒内の流れに影響を与えたからだと考えられる．

3-(a) 0D
3-(b) 4D
Fig.3. Visualization result (Velocity magnitude contour plot)
モニタ点から得たエンジン筒内の流速の時系列データを図 4 に示す．各グラフは横軸がクランク角度(CA)を表しており，圧縮上死点時を 0 deg としている．
図 4 は 0D モデル，4D モデル共に約 50 サイクル分のアンサンブル平均値も示している．圧縮行程で流速の平均値からのばらつきがあることがわかった．
0D モデルと 4D モデルの流速データを用いてサイクル変動成分を抽出した．サイ
クル変動成分の抽出は次の式を用いた[1]．𝑢𝑐 = |𝑢̅ − 〈𝑢〉| (𝑢𝑐 :サイクル変動成分，𝑢̅:時間平均流速，〈𝑢〉:アンサンブル平均)．ここでの時間平均流速は各サイクルの流速の時
系列データに LPF を掛けたものとした．
図 5 に代表サイクル(10，20，30，40，50 サイクル目)のサイクル変動成分の比較を示す．サイクル変動成分は 0D モデル，4D モデルともに吸気行程時が最も大きくなっており，4D モデルは顕著に大きくなっている．これは吸気管に取り付けた乱流促進体の影響だと考えられる．圧縮工程においてサイクルごとに流速のばらつきが見られたが，サイクル変動成分は小さくなっていた．これより，吸気工程ではサイクル変動成分が多く含まれており，圧縮工程では乱流成分が多くなることがわかった．

102

JSST 2022

4-(a) 0D

4-(b) 4D
Fig.4. Comparison of Velocity magnitude

5-(a) 0D (Inside cylinder)
5-(b) 4D (Inside cylinder)
Fig.5. Comparison of cycle-to-cycle variation

4. まとめ

延長吸気管付きのモデルを数値計算することができた．また，延長吸気管がついていないモデルとの比較を行い，吸気面での乱流分布がエンジン筒内の流れ場に影響を与えることがわかった．
得られた時系列データより，サイクル変動成分を抽出することができた．これにより，吸気面での乱流分布が影響を与えていることがわかった．

参考文献[1] Satoshi HOKIMOTO, Tatsuya KUBOYAMA, Yasuo MORIYOSHI, Takahiro
WATANABE, Minoru IIDA: YAMAHA MOTOR TECHNICAL REVIEW, Analyses ofcycle-to-cycle variation of combustion in low load conditions by visualization of fuelmixture formation processes using an optical engine, (2016), 95-101.
[2] 溝渕泰寛: 自動車技術会 2019 春季大会フォーラム, エンジンのモデルベース開発を加速する国産燃焼ソフトウェア HINOCA, 20194415(2019), 5-9.

103