JSST2021
The 40th JSST Annual International Conferenceon Simulation Technology

Student Session Proceedings

Doshisya University
Kyoto, JAPAN
September 1–3, 2021

Supporting Companies

iiSession ScheduleDay 1 (Wednesday, 1st September, 2021)
Time

Hall A

Hall B

Hall CHall D

09:30 – 09:45

Opening Ceremony (Hall A)

09:45 – 10:45

Plenary Talk 1 (Symposium 1, Hall A)

10:45 – 11:00

Group Photo (Hall A)

11:00 – 12:15

Symposium 1 1
Invited Talk 1
Tutorial Talk 1
Session 1 (2)Hall E

OS8 1

OS3 + OS9 1

OS2 1

OS7 1

Session 5 (4)Session 8 (3)Session 12 (3)Session 16 (3)Lunch Time

12:15 – 13:30

Corporate Seminar 1 (Hall A)
13:30 – 14:15

Differential equations and machine learning/neural networks
Koji Maruyama (Wolfram Research)

14:30 – 15:30

15:45 – 16:45

17:00 – 18:00

Symposium 1 2

OS8 2

OS3 + OS9 2

OS2 2

OS7 2

Session 2 (3)Session 6 (3)Session 9 (3)Session 13 (3)Session 17 (3)Symposium 1 3

OS8 3

OS3 + OS9 3

OS6 1

OS11 1

Session 3 (3)Session 7 (3)Session 10 (3)Session 14 (3)Session 18 (3)Symposium 1 4

Symposium 2 1

OS6 2

OS11 2

Session 4 (3)Session 11 (3)Session 15 (2)Session 19 (2)iiiDay 2 (Thursday, 2nd September, 2021)
Time

Hall A

11:30 – 12:30

12:30 – 13:30

Hall CHall DHall E

Plenary Talk 2 (Symposium 2, Hall A)

09:00 – 10:00

10:15 – 11:15

Hall B

Symposium 1 5
Session 20 (3)Symposium 2 2
Invited Talk 2
Invited Talk 3
Session 22 (2)Symposium 3 1

OS11 3

OS1 1

Session 24 (3)Session 26 (3)Session 28 (3)Symposium 1 6

Symposium 2 3

Symposium 3 2

OS10 1

OS1 2

Session 21 (3)Session 23 (2)Session 25 (3)Session 27 (3)Session 29 (3)Lunch Time
Corporate Seminar 2 (Hall A)

13:30 – 14:15

The Trends of Latest Personal Supercomputer: Performance and App Support Status
Keiou Nishimura and Rika Yamashita (Visual Technology)

14:30 – 15:00

Plenary Talk 3 (Symposium 3, Hall A)

15:00 – 15:05

Announcement about JSST2022 (Hall A)

15:15 – 16:15

Shotgun Presentations of Student Session (Hall A)

16:30 – 18:00

Student Session (Poster Session, Remo)iv

Day 3 (Friday, 3rd September, 2021)
Time

Hall A

11:30 – 12:45

Symposium 3 3

OS4 1

RS 1

Session 30 (3)Session 33 (3)

OS4 2

RS 2

Session 31 (2)Session 34 (2)Hall DHall E

OS10 2

OS5 1

Session 39 (3)Session 42 (3)Symposium 3 4

OS10 3OS5 2

Session 37 (3)Session 40 (3)Session 43 (3)Tutorial Talk 2
Session 36 (2)Lunch Time

12:45 – 13:30

13:30 – 14:15

Hall CPlenary Talk 4 + Invited Talk 4 (Symposium 3, Hall A)

09:00 – 10:00

10:15 – 11:15

Hall BOS4 3

Symposium 2 4

Symposium 3 5

OS10 4OS5 3

Session 32 (2)Session 35 (2)Session 38 (2)Session 41 (2)Session 44 (2)14:15 – 14:30

Break

14:30 – 15:00

Closing Ceremony (Hall A)v

Symposium
Plenary Talk
Presenter:

Prof. Hiroaki Natsukawa
Academic Center for Computing and Media Studies,
Kyoto Universit

Title:
Understanding System Dynamics by Combining Data-Driven Analysis andInformation Visualizatio

Abstract:
Natural systems are often complex and dynamic. Reductionism is notuniversally applicable for natural complex systems found in biology andelsewhere where behavior is driven by complex interactions between many
components acting together in time nonlinear dynamic systems. In this talk, I
will introduce a minimalist paradigm, empirical dynamic modeling (EDM), forstudying non-linear systems. It is a data-driven approach that uses timeseries information to study a system by reconstructing its attractor (ageometric object) that embodies the rules of a full set of equations for thesystem. Particularly, I will introduce a few recent works of visual analytics tosupport the identi cation and mechanistic interpretation of system statesusing integration EDM with visualization techniques, leads to ourunderstanding of complex dynamics

Biography:
Hiroaki Natsukawa is a junior associate professor / senior lecturer in
Academic Center for Computing and Media Studies, Kyoto University. He
received a Ph.D. in engineering from Kyoto University in 2013. Althoughformally trained as a researcher in the eld of biomedical engineering at
Kyoto University, He has successfully crossed elds into other areas such as
information visualization, neuroscience, and more recently into nonlinear
dynamics. Currently, he has worked in the eld of information visualizationand his work focuses on developing visual analytics systems enabling datavi

driven analytical reasoning by empirical dynamic modeling in collaboration
with UC San Diego. He is currently an associate editor of Journal of
Visualization, was a The 36th Annual Meeting of Japan Biomagnetism andBioelectromagneics Society Local Committee, IEEE Paci cVis 2019 Poster
Co-chairs, IEEE Paci cVis 2018 Publication Chair, The 47th Symposium of
Visualization Society of Japan Organizer, NICOGRAPH International 2017
Local Committee Chair, JSST2016 Publication Chair, JHBM18 Local
Committee.

Invited Talk
Presenter:

Prof. Tomohiro Sogabe
Department of Applied Physics, Nagoya Universit

Title:
Computation of Singular Values for Generalized Tensor Su(Joint work with Asuka Ohashi (National Institute of Technology, Kagawa

Abstract:
We consider computing singular values of the generalized tensor sum of thefor

T := In ⊗ Im ⊗ A + In ⊗ B ⊗ Iℓ + C ⊗ Im ⊗ Iℓ,(1

where is the × identity matrix and A ∈ R ℓ×ℓ, B ∈ R m×m, A ∈ R n×n . The
mathematical symbol ⊗ denotes tensor product (or Kronecker’s product). A
simple example of the tensor product is given below.

a11b11a11B a12 B
a21b11A ⊗ B :==
[a21B a22 B]
a11b21
a21b21vii

a12b11
a22b11
a12b21
a22b21

a11b12
a21b12
a11b22
a21b22

a12b12
a22b12
a12b22
a22b22

where A := a
. Similarly, A ⊗ B ⊗ C is calculated recursively by[ 21 a22]

a11 a12(A ⊗ B) ⊗ C or (A ⊗ B) ⊗ C . The size of the generalized tensor sum in (1)is of n 3 × n 3, which can be extremely large even if matrices , , are small.
Indeed, for , , being 1,000 × 1,000 matrices, the generalized tensor su
can be a 1,000,000,000 × 1,000,000,000 matrix
Though the generalized tensor sum can be very large, it is easy to compute
some fundamental quantities of linear algebra, such as determinant andeigenvalues. In fact, the eigenvalues ofcan be written as the sum of
eigenvalues of relatively much smaller matrices , , andthan . On theother hand, unlike eigenvalues, it is dif cult to compute the singular values osince there is no such a simple relation between the singular values oand the singular values of , , andIn this talk, we present recent progress on ef cient computational algorithms[1,2,3] for some speci c singular values of the generalized tensor sum ,which are based on the use of the notion of numerical multilinear algebra
Reference[1] A. Ohashi, T. Sogabe, On computing maximum/minimum singular valuesof a generalized tensor sum, Electronic Transaction on Numerical Analysis,
43 (2015), pp. 244-254[2] A. Ohashi, T. Sogabe, On computing the minimum singular value of atensor sum, Special Matrices, 7 (2019), pp. 95-106[3] A. Ohashi, T. Sogabe, On computing an arbitrary singular value of a tensor
sum, in preparation

Biography:
Tomohiro Sogabe was graduated from the department of Applied Physics, theUniversity of Tokyo in 2001 and received Ph.D. from the same university in2006. He worked at Nagoya university as an assistant professor, and Aichi
prefectural university as an associate professor. He is currently an associate
professor at the department of Applied Physics, Nagoya University, Japan. He
served as editors of JSIAM Letters, the Transactions of JSIAM, and currently
he serves as an editor of JJIAM journal, Springer. He will be a member of theboard of directors of JSIAM from June 2021. His research interests include
numerical linear algebra, numerical multilinear algebra, and scienti ccomputing. He published more than 70 international journal papers and 15
domestic journal papers include applied mathematics journals, computational
physics journals and quantum computing journal, such as AMC, AML, ETNA,
JCAM, JJIAM, NLAA for applied mathematics; IEEE TMTT, JCP, PRB, PRE
for computational physics; QIC for quantum computing journalviii

Tutorial Talk
Presenter:

Prof. Kuniyoshi Abe
Faculty of Economics and Information, Gifu Shotoku
Gakuen Universit

Title:A Numerical Study of Parallel Variants of GPBiCG
Method with Stabilization Strategy for Solving Linear Equation(Joint work with Soichiro Ikuno and Gerard L.G. Sleijpen

Abstract:
The hybrid Bi-conjugate gradient (Bi-CG) methods such as Bi-CG stabilized(Bi-CGSTAB), Generalized Product-type based on Bi-CG (GPBiCG) are wellknown for ef ciently solving linear equations, but we have seen theconvergence behavior with a long stagnation phase. In such cases, it isimportant to have Bi-CG coef cients that are as accurate as possible. We
introduce the stabilization strategy for improving the accuracy of the Bi-CG
coef cients.
On present petascale high-performance computing hardware, the main
bottleneck for ef cient parallelization is the inner products which require a
global reduction. The parallel variants of Bi-CGSTAB reducing the number ofglobal communication phases and hiding the communication latency havebeen proposed. In this paper, therefore, following the analogy of Cools et al.,we design parallel variants of GPBiCG, and examine the convergence speedof the parallel variants with the stabilization strategy

Biography:
Kuniyoshi Abe obtained Ph. D at Nagoya University in 1999. He began his
career as Assistant professor at Anna National College of Technology in
1998. He moved to Advanced Computing Center, RIKEN as a Contract
researcher in 1999. He was associate professor at Faculty of Economics andInformation, Gifu Shotoku University in 2003, and is professor from 2012. His
research eld is numerical linear algebra, especially, he is interested in fastand iterative solvers for linear equationsix

Symposium2
Plenary Talk
Presenter:

Prof. Hajime Ishihara
Department of Materials Engineering Science, Osaka Universit

Title:
Interaction of Optical Vortex and Matter Systems: from Microscopic to
Macroscopic Regimes and Linkage between The

Abstract:
A remarkable effect in the interaction between the optical vortex and mattersystems is the transfer of the orbital angular momentum. A lot of theoreticaland experimental demonstrations have been conducted for macroscopic
rotational manipulation of matters from microscale to nanoscale objects. On
the other hand, in the microscopic regime, the optical vortex changes theselection rule of electronic transitions, which induces unconventional spin
dynamics in solids. Interestingly, the macroscopic and microscopic degrees of
freedom are sometimes linked to each other in a unique way via linear andnonlinear optical responses. In this talk, I will present some examples of theabove topics and puzzling aspects about them

Biography:
Hajime Ishihara is a Professor in Graduate School of Engineering Science at
Osaka University. Since 2016, he has been the project representative of
JSPS Grand-in Aid for Scienti c Research on Innovative Area “Nano-Material
Manipulation and Structural Order Control with Optical Forces”. His researchinterests are the microscopic interaction between light and nanostructures
and the optical functionalities arising from their interaction. He is the fellow of
the Japan Society of Applied PhysicsxInvited Talk
Presenter:

Prof. Yoshiki Kohmura
RIKEN SPring-8 Cente

Title:
Generation of X-ray Vortices by Bragg Re ection from Crystals andMicroscopy to Diagnose Topological Charge Distributio

Abstract:
We developed a method to transfer the structure of dislocation in a crystal to
the Bragg X-ray wave front and performed a phase measurement that could
diagnose the dislocation. This method helps to enhance the functionality ofmaterials affected by dislocations. The two-beam X-ray interferometry was
performed to measure the Bragg re ected X-ray wave front from a spiraldislocation area in a SiC crystal. As a result, a pair of fork patterns wasobserved. The number of fringes changed by +1 and -1 at the phaseanomalies due to generation of an X-ray vortex with the topological charge of+1. Kinematic diffraction simulation showed that X-ray vortices are generated
on Bragg re ected X-rays from spiral dislocation areas. This result is realized
by using suf ciently large spatial coherence of the irradiated X-rays

Biography:
Yoshiki Kohmura is currently a Team Leader of Synchrotron Radiation
Imaging Instrumentation Team at SPring-8 facility, Japan. He received his
PhD in the Department of Physics, Graduate School of Science, theUniversity of Tokyo in 1994. He has been working for RIKEN SPring-8 Center
since 1996. He is involved in developing various x-ray optical elements andmethodologies for synchrotron radiation experiments. In the past ten years,
he studied the novel waveguiding phenomenon of x-rays occurring insideslightly deformed crystals and its application to fast switching of x-raysxi

Invited Talk
Presenter:

Dr. Toru Tsujimura
National Institute for Fusion Science, National
Institutes of Natural Science

Title:
Propagation Properties of Optical Vortex in Magnetized Plasm

Abstract:
We investigated theoretically propagation properties of an optical vortex witha helical wavefront in cold uniform magnetized plasmas in an electron
cyclotron (EC) range of frequencies. We describe the effects of the helical
wavefront of the vortex EC wave on the wave elds in magnetized plasmas
as anisotropic media. These effects become signi cant as the topologicalcharge of the vortex EC wave increases or the distance from the phasesingularity becomes small. The different propagation properties are also
con rmed in the propagation of Laguerre–Gaussian beams on threedimensional simulations with the nite element method. We found aremarkable phenomenon on mode conversion in the vortex EC wave
propagation, which suggests an optical vortex can be a new tool to heat highdensity plasmas

Biography:
Toru Ii Tsujimura received the B.Eng., M.Eng., and D.Eng. Degrees from the
University of Tokyo in 2009, 2011, and 2014, respectively. Since 2014, he has
been working as an assistant professor at National Institute for Fusion
Science, National Institutes of Natural Sciences (NINS). His current research
focuses on wave heating in fusion plasma at the Large Helical Device. He
was a representative of NINS programs for cross-disciplinary study in 2018and 2019, leading the program on helical wavefront measurements of opticalvortices emitted by high-energy electronsxiiSymposium3
Plenary Talk
Presenter:

Prof. Kenichi Yoshikawa
Self-organization Science Research Center, Doshisha
University / Institute for Advances Study, Kyoto Universit

Title:
Cooperation of Real-World Modeling with Simulation toward the Problem“What is Life?

Abstract:
Living matter is manipulating large number of different molecules to create
spatio-temporal self-organization. Currently, many numerical studies havebeen carried out to shed light on the underlying mechanisms of life.
Sometimes, these studies remain as a possible explanatory interpretation
because of the complexity of actual living systems. In the present talk, I would
like to show some examples of real-world modeling, which help us to obtaindeeper understanding in biological systems. 1) Emergence of cell-likestructure and function under crowding condition with macromolecules. 2) Selforganization of cellular assembly. 3) Self-generation of macroscopic regularmotion for the assembly of nano-scaled self-propelling objects

Biography:
Name
Kenichi Yoshikawa
URL: http://dmpl.doshisha.ac.jp/
Academic Degree
PhD, Physical Chemistry, Graduate School of Engineering, Kyoto University(1976)

xiii

Professional Appointments
1998 – 2012
Professor, Department of Physics, Graduate School of Science, Kyoto
University
2011 – 2014
Chair, Commission of Biological Physics C6, IUPAP
2012 –
Professor, Faculty of Life and Medical Sciences, Doshisha University
2019Specially Appointed Prof., Institute for Advances Study, Kyoto Universit

Plenary Talk
Presenter:

Prof. Helmut Schiessel
Cluster of Excellence Physics of Life, Technical
University Dresde

Title:The mechanical genom

Abstract:
In addition to the classical layer of genetic information, DNA molecules also
carry other layers of information. Two possible layers are discussed. One is amechanical layer, namely that the base pair sequence in uences themechanical and geometrical properties of the DNA molecules which in turnmight guide their packing inside the cell. A second possible layer is the speedof translation in the protein-producing ribosomes where the speed in uences
the quality of the protein product. We demonstrate that the degeneracxivBiography:
Helmut Schiessel studied physics in Freiburg where he did his Ph.D. intheoretical polymer physics in 1997. He then worked as a postdoc at two
Universities of California, in Santa Barbara and in Los Angeles. In 2000 hejoined the MPI for Polymer Research in Mainz where he was in charge of abiophysics research project. From 2005 to 2020 Schiessel headed the chairof Theoretical Physics of Life Processes at the Lorentz Institute, Leiden
University. Since January he heads the group Theoretical Physics of Living
Matter at the Cluster of Excellence Physics of Life.

Invited Talk
Presenter:

Prof. Yoshiteru Yonetani
National Institute for Quantum and Radiological
Science and Technolog

Title:
Water on the DNA Surface: Microscopic Insight from Molecular Dynamics
Simulation

Abstract:
Water on the DNA surface exhibits characteristic behavior such as ice-likewater network along the narrow groove of DNA. By using molecular dynamics
simulations, we obtained detailed picture for such water; our characterization
fully describes the water network patterns and their sequence variation[Biophys. J. 97, 1138 (2009)]. Analysis on such water network further
revealed microscopic origin of the slow water-exchange kinetics [Biophys.
Chem. 160, 54 (2012)]. As another relevant topic, I will also mention a relation
with DNA damage production. Our recent study [Chem. Phys. Lett. 749,
137441 (2020)] found that site-dependent water accessibility for DNA
backbone is well correlated with the probabilities of DNA damage production
by OH radical attacking. This suggests that the site-dependent hydration
property is a key factor for DNA damage production

xv

Biography:
Yoshiteru Yonetani received B.E. and Dr. degrees from Keio University in1998 and 2002, respectively. He worked as a postdoctoral researcher at Nara
Women’s University (2002-2003) and Japan Atomic Energy Agency(2004-2008). He is currently a senior researcher in National Institutes of
Quantum and Radiological Science and Technology. His current interests
include molecular dynamics and quantum dynamics in biological andchemical physics

Tutorial Talk
Presenter:

Dr. Miyuki Sasaki
Collaborative Laboratories for Advanced
Decommissioning Science, Japan Atomic Energy
Agenc

Title:
Visualization of Dose Rate Distribution around Fukushima Daiichi Nuclear
Power Plant using Arti cial Neural Networks

Abstract:
This study proposes a method of visualizing the ambient dose rate
distribution using arti cial neural networks (ANN) from airborne radiation
monitoring results. The ANN method was applied to the results of the airborne
radiation monitoring which was conducted around the Fukushima Daiichi
Nuclear Power Plant by an unmanned aerial vehicle. The ANN wasconstructed by training data consisting of input variable dataset (radiation
count rate, altitude, topographic data, photographic RGB data) and objective
variable dataset (the air dose rate data at 1 m above the ground level). The
reliability of the ANN method was evaluated by comparison with the groundbased survey data. The dose rate map created by the ANN method
reproduced ground-based survey results better than traditional methods

xviBiography:
Miyuki Sasaki received Dr. Eng. degrees from the University of Nagoya,Japan, in 2021. She became a researcher at the Japan Atomic Energy
Agency in 2015. Her current research interests include radiation
measurementxviiCommittee
General Chair:
Hiroaki Nakamura (National Institute for Fusion Science
Conference Chair:
Takahiro Kenmotsu (Doshisha University
Conference Co-Chair:
Susumu Fujiwara (Kyoto Institute of Technology)
Katsuhiko Yamaguchi (Fukushima University
Program Chair:
Seiki Saito (Yamagata University)
Hiroshi Tamura (Chuo University
Program Co-Chair:
Nobuaki Ohno (University of Hyogo)
Yoshihide Shibata (National Institute of Technology, Gifu College )Yoshihisa Fujita (Ritsumeikan University
Post-Conference Chair:
Susumu Fujiwara (Kyoto Institute of Technology
Post-Conference Co-Chair:
Tomoko Mizuguchi (Kyoto Institute of Technology)
Hiroaki Ohtani (National Institute for Fusion Science)
Daisuke Ishihara (Kyushu Institute of Technology
Publication Chair:
Taku Itoh (Nihon University)
Publication Co-Chair:
Hiroto Tadano (University of Tsukuba)
Soichiro Ikuno (Tokyo University of Technology)
Finance Chair:
Hiroki Sakuta (Doshisha University
Publicity Chair:
Satoshi Takatori (Doshisha University
Exhibition Chair:
Yoji Matsumoto (Cybernet Systems Co., Ltd.
Conference Adviser:
Amane Takei (University of Miyazaki)
Takuo Yasunaga (Kyushu Institute of Technology)
Motoi Wada (Doshisha University)

xviii

Table of Contents for Student Session
Day 2 (Thursday, 2nd September, 2021)
Student Session (Poster Session)
P01 Defect Detection by Simulation and Machine Learning .........................................1
Yu Shirafuji, Tomonori Yamada
P02 Short-term prediction using RBF network with coordinate transformation .................. 5
Kohei Saito, Satoshi Kitayama
P03 Discrimination of Japanese drum player using machine learning ..............................9
Kanako Toyoda, Hirosh Tamura
P04 Perception of Direction of Moving Sounds. -Study on Adaptive Method- ................. 13
Tetsu Takizawa, Akihiro Kudo
P05 A Method for setting up Servers for Edge Computing in Fifth-Generation Mobile
Communication Systems ........................................................................... 17
Ryota Ehara, Hiroshi Tamura
P06 Solving Sudoku with Graph Convolutional Networks ......................................... 21
Koichiro Ishii, Hiroshi Tamura
P07 On delivery methods using UAV in hilly and moun- tainous areas ...........................25
Yuto Mashiko, Hiroshi Tamura
P08 Development of visualization system for radiation by mixed reality technolog ............29
Masato Iida, Sotaro Kudo, Rika Kawahara, Katsuhiko Yamaguchi, Seiki Saito
P09 Computation of Micro-particles and Fluid Flow through Porous Structure .................32
Nobuyuki Hirooka, and Satoru Ushijima
P10 UAV assisted Survivor Search Scheme with Survivor Density Distribution Prediction ...36
Chen Haosong, Kazutoshi Yoshii, Shigeru Shimamoto
P11 Elastodynamic analysis of elastic wave scattering using coupling method of CQBEM andFEM ................................................................................................... 40
Haruhiko Takeda, Takahiro Saitoh
P12 Particle-scale FSI computation for internal ﬂuidization in gravel-particle bed by upward
water jet ............................................................................................... 42
S. Ushijima, J. Ohno, D. Toriu, Y. Ueno
P13 Embodiment Design of Biped Robot for Walking on Rough Ground ....................... 46
Kakeru Yokoyama, Van-Tinh Nguyen, Hiroshi Hasegawa

xix

P14 Effect of Adsorption and Reactive Behavior of Water Vapor in Reactive Magnetron
Sputtering Simulations using Berg’s Model ..................................................... 49
Allen Vincent B. Catapang, Motoi Wada
P15 Target ionization and surface damage models of graphite via nanosecond pulse
irradiation ............................................................................................. 53
James Edward Hernandez, Kengo Moribayashi, James Koga, Motoi Wada
P16 Analysis of the space charge effect in an ion mobility spectrometer .........................57
Keith Nealson M. Penado, James Edward A. Hernandez, Allen Vincent B. Catapang, Motoi
Wada
P17 Molecular Dynamics Simulation of Clustered DNA Damage Composed of
Apurinic/Apyrimidinic Sites ....................................................................... 61
Kazushi Terakawa, Kawanami Ryuta, Haolun Li, Susumu Fujiwara, Tomoko Mizuguchi,
Yoshiteru Yonetani, Hiroaki Nakamura
P18 Consideration of environment-friendly manufacturing in China ............................. 65
Zhou Peihan, Hitoshi Sugimoto, Akira Otsuka
P19 The study on effects of wireless standards on dense Wi-Fi environments with NS3 ...... 69
Hiroshi Iwata, Sakie Horiuchi
P20 Computational method for interactions between deformable objects and ﬂuid ﬂows using
immersed boundary method and mass spring model ........................................... 73
Niku Guinea, Daisuke Toriu, Satoru Ushijima
P21 Simulation of elastic waves in micropolar bimaterials using 2-D M-EFIT ................. 77
Yusuke Suzuki, Takahiro Saitoh
P25 Molecular Dynamics Study of Water Dynamics in Zwitterionic Polymer Brush-Water Interface ..................................................................................................79
Yuya Fujinaga, Susumu Fujiwara, Tomoko Mizuguchi
P26 Development for High-Speed of High-Frequency Electromagnetic Fields Solver: ADVENTURE FullWave ..................................................................................... 81
Kento Ohnaka, Amane Takei
P27 Development and application for wave-sound ﬁeld Solver: ADVENTURE Sound .......85
Ryo Yoshidome, Amane Takei
P28 Molecular dynamics simulation and numeric calculation of a damaged polyethylene assuming tritium substitution and decay ................................................................. 88
Ryuta Kawanami, Susumu Fujiwara, Hiroaki Nakamura, Kazumi Omata
P29 Magnetization process analysis of magnetic material with void defect ......................92
Yasuyuki Igarashi, Katsuhiko Yamaguchixx

P30 Mechanical properties of hydrogen embrittled austenitic stainless steel .................... 95
Keisuke Uchida, Katsuhiko Yamaguchi
P31 Simulation for estimation of radiation source distribution measured by NaI(Tl) scintillation
detector ................................................................................................ 99
Sota Suga, Haruka Yamazaki, Katsuhiko Yamaguchi
P32 Simulation of thermoluminescence measurement of structural materials ..................102
Misaki Uesugi, Tsugiko Takase, Katsuhiko Yamaguchi
P33 Study for Development and application of Nonsteady High-Frequency Electromagnetic
Field method ........................................................................................ 105
Souta Nishimura, Amane Takei
P34 Localization phenomenon of large and small particles by congestion degree and boundary
conditions : Cell model experiment on cm scale .............................................. 108
Maho Kuroda, Satoshi Takatori, Chwen-Yang Shew, Takahiro Kenmotsu,
Kenichi Yoshikawa
P35 Sightseeing Guidance System to Maximize Satisfaction Using Real-Time Spot
Information ..........................................................................................112
Yuya Sato, Yoko Nakajima, Hirotoshi Honma
P36 Evaluation of solvent-accessible surface area of back- bone hydrogen in telomeric DNA for
studying tritium resistance of DNA ..............................................................114
Yohei Tsuchida, Seiki Saito, Hiroaki Nakamura, Yoshiteru Yonetani, Susumu Fujiwara
P37 Smart Diagnosis on Disease State based on Characteristic Cracking Pattern by Stretching
Response of Tissue Slice .......................................................................... 117
Tagami Yukiho, Satoshi Takatori, Takahiro Kenmotsu, Kenichi Yoshikawa

xxi

JSST 2021

シミュレーションと機械学習による内部欠陥判定Defect Detection by
Simulation and Machine Learning
白藤優 1*,山田知典 1
Yu Shirafuji1*,Tomonori Yamada1

1
1

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

Department of System Innovation, Graduate School of Engineering, University of Tokyoshirafuji-yu@g.ecc.u-tokyo.ac.jp

Abstract. It is very important to detect defects in structures to prevent accidents related tothem. A conventinal way depends on experts’ experiences and skills. Therefore some researches where people try to create a system to detect defects by using machine learning havebeen conducted. However, it is diffuclt to collect sufficient data points which you use formachine learning in the real world. In our research, I created data points by using computer
simulation and completed machine learning by using these data points. Furthermore, considering cases where simulations require huge computation costs, we conducted numerical experiments to find out a way to create highly accurate machine learning methodology withsmall data points. Concretely we introduced the theory of design of experiments to a sampling of machine learning and confirmed an improvement of the accuracy of the system.
Keywords: Computer simulation, Machine learning, Design of experiments

1. はじめに様々な構造物の事故を防ぐには構造が有する欠陥を未然に発見することが重要になる 1)。従来、このような欠陥発見のために打音検査が用いられることが一般的であるが、これは専門家の経験に頼った方法であるため、近年このような非破壊検査に機械学習の技術が利用されている 2-3)。しかし機械学習による欠陥判定システムには十分な学習データが必要になる一方で欠陥を有する構造についてデータを収集することは難しい。そのためシミュレーションによって学習データを作成することでこの課題を解決する研究が行われている。つまりシミュレーションでデータセットを作成し、それらを学習に用いることで、観測することができる情報から構造物の内部欠陥情報の推定システムを作成することを目標とした研究である。しかしシミュレーションに大きな計算量と時間がかかることはこのアプローチの課題だ。そのため可能な限りシ

1

JSST 2021
ミュレーションにより作成するデータ数を減らすことが求められている。本研究では、
シミュレーションと機械学習を利用して欠陥情報を推定するシステムを作成し、そのシステムを利用して実験を行うことで少ないデータ数で高精度のシステムを作成す
る方法について考える。また行うシミュレーションについてはまずは予備的な検討として、2 次元の定常熱伝導の温度分布シミュレーションとした。

2. 内部欠陥を有する構造のモデル化
本研究では単純化のため構造、構造が有する欠陥どちらについても正方形とした。
また本研究では構造に対して定常熱伝導シミュレーションを行うため、メッシュを作成する必要がある。そのため正方形上に 2500 個の点を取りこれらの点を頂点とする
直角三角形を１要素として切り分けた。またこれらの点を頂点として、大きさ 1 から
大きさ 5 の正方形欠陥を作成した。そのため図 1 のように内部に 1m 四方から 5m 四方の正方形欠陥を有する 49m 四方の正方形構造を作成した。

3. 定常熱伝導シミュレーション
本研究では有限要素法による定常熱伝導のシミュレーションを行うことで温度分布データを作成するプログラムを実装した 4)。境界条件は左辺の中心を 100℃、右辺を 0℃とした。これを 2 章で作成した構造に対して実行することで、各構造の温度分
布データを作成した。またデータのまとめ方について、入力となる温度分布について
は今回観測できる表面と考えられる左辺一辺に注目し、欠陥がない場合の左辺一辺との温度差をまとめた。また出力については 1m 四方から 5m 四方の欠陥の大きさと位置を示す欠陥中心の座標をデータセットにまとめた。

図1

定常熱伝導シミュレーション実行例

4. 機械学習本研究では特徴量の抽出に優れた畳み込みニューラルネットワーク(CNN)のモデルを実装した 5,6)。このモデルを利用して欠陥の位置と大きさに関する推定を回帰問題として取り組んだ。本論文では大きさ推定の説明に集中する。またデータについては
3 章で作成したデータセットを 4:1 に学習データと訓練データにランダムで分け学習
を行った。

2

JSST 2021
4.1 5 層ニューラルネットワークモデルによる学習結果

まずはモデルを、入力層、畳み込み層、プーリング層、全結合層、出力層からなる
単純なものとし、学習を行った。正答とする基準を誤差 0.5 以下、誤差が欠陥の大きさの 20%以下とした場合の正答率を算出したところ、基準を絶対値とした場合、欠陥が大きいほど精度が落ちていること、基準を相対値とした場合、欠陥が小さいほど精度が落ちていることが確認できた。
4.2 9 層ニューラルネットワークモデルによる学習結果
層数にルールはないが一般的に層を増やせば精度が向上することから 7)畳み込み層とプーリング層を 3 セット重ね、モデルを深くし同様に学習を行った。その時の結果を
先程の結果と比較したものが以下の図 2,図 3 である。正答とする基準がどちらの場
合も 4.1 のモデルよりも十分精度が上がっていることが確認できる。

図 2 大きさ毎正答率比較：誤差 0.5 以下 図 3 大きさ毎正答率比較：誤差 20%以下
4.3 データを減らした場合
これまでは全てのデータを用いて学習を行ってきたが現実的ではないため、データ
数を 250 から 5000 まで 20 段階で変更しその時 loss を比較した。またこの際、実行時間削減のためはじめに作成した単純なモデルを利用した。その時の loss の経過が下の図 4、最終的な loss の値が下の図 5 のようになった。これらからデータ数が少なくなると誤差が大きくなる傾向があること、またデータ数 250 では大きく精度が落ちていることが確認できる。

図 4 各データ数の loss 経過：大きさ推定 図 5 各データ数の最終的な loss：大きさ推定

5. 実験、考察
前章で利用するデータを少なくすると精度が落ちること、またデータ数が 250 の時は特に精度が落ちていることが確認できた。本章では少ないデータ数でも精度を保つための方法を見つけるため実験を行う。具体的には学習に用いるデータをランダムではなく、実験計画法に基づいて選択することによって効率的に学習を進め、結果として高精度の推定システムが作成できるのか確認する。
5.1 実験方法
前述の通りこの実験ではデータ選択に実験計画法を用いることで少ないデータ数で高精度が実現できるのか確認する。実験計画法とは効率の良い実験方法を設計し結果を適切に解析することを目的としており、各変数が結果にどのように影響しているかという関係を効率的に求めるための手法である。今回はこの実験計画法の直行配列

3

JSST 2021

実験という手法を利用する。これは実験の組み合わせを工夫することで少ない実験数
で各因子の結果への影響を知る手法である。今回の問題は因子については欠陥の大きさ、欠陥の x 座標、y 座標の３つ、水準については欠陥の大きさ 1m 四方から 5m 四方、位置についても 5 分割した時のどこに属するかと考えることで 5 つと考えることできるため、3 因子 5 水準の直行配列実験を応用した。つまり 5 水準直交表に基づいて学
習データを選ぶことで精度の向上を目指す。その他の設定については 4.1 で単純なモデルで学習した際と同様にし、その場合の結果と比較する。
5.2 実験結果
実験計画法に基づいて 250 から 1000 個のデータを選び学習を行った際の loss と4.1 の loss を比較したものは下の図 6 のようになっている。この結果から少ないデータ数で loss が小さくなっていること、またデータ数が 250 の際には大きく loss が小さくなっていることが確認できる。

図6

実験計画法 loss 比較：大きさ推定

6. おわりに
シミュレーションにより作成したデータセットを機械学習に用いることで構造の内部欠陥判定システムを作成した。また本システムを利用した実験から実験計画法に基づいたデータ選択が少ないデータでの効率的な機械学習に有効であることを示した。

参考文献[1] 田中正吾,山田実『信号伝播モデルに基づく電磁波レーダによるコンクリート構
造物の非破壊検査』,計測自動制御学会,39(5),pp.432-440[2] 高橋修司,宮島雅弥,堀口敦史,中上京治,茂木和弘,白石洋一,須田高史『機械学習を用いた打音による鋼管柱の非破壊欠陥推定』,検査技術,23(10).pp.15-23[3] 堀辺忠志,岡村弘之『ニューラルネットワークによるき裂を有するはりの逆問題
解析』,日本機械学会,64(617).pp.1-6[4] 邵⻑城『基本からわかる有限要素法』森北出版株式会社,2008[5] Keras,https://keras.io/,2021.07.01.
[6] Diederik Kingma; Jimmy Ba Adam: A Method for Stochastic Optimization,
arXiv:1412.6980v1,2021.07.01[7] 岡谷貴之『画像認識のための深層学習の研究動向』,人工知能学会,31(2),pp.169-179

4

JSST2021

Short-term prediction using RBF networkwith coordinate transformation
Kohei Saito1*, Satoshi Kitayama2
1

Graduate school of natural science and technology, Kanazawa University
2
Advanced manufacturing technology institute, Kanazawa University
ksaitou.cl.816@stu.kanazawa-u.ac.jp

Abstract. In this paper, we present a short-term prediction using RBF network with affine
transformations. By incorporating an affine coordinate transformation into the conventional
RBF network, we can make predictions that take into account the direction of the data. In otherwords, the prediction is performed in a rotational coordinate system and the predicted value are
returned to an absolute coordinate system. It is found from numerical result that the prediction
accuracy is improved by using affine transformation.
Keywords: Short-term prediction, RBFnetwork, Time series, Affine transformation

1. Introduction
RBF network is one of the artificial neural networks that use Gaussian functions as basis
functions. Short-term prediction by RBF network has been shown to be highly accurate andcomputationally inexpensive(1)(2)(3). However, short-term prediction using RBF networks does
not take into account the direction of the time series. So, we introduce a method considering
the direction of the time series by affine transformation.

2. Short-term prediction of time series using RBF network
Nonlinear modeling with RBF networks 𝑦" (") can be expressed as follows(

𝑦" (") = $

("%&)

$)*

𝑤$

("%&)

𝐾 '𝒙("%&) , 𝒙𝒋

(1)*

Where 𝑡 is the current time, 𝜏 is the delay time, 𝑚 is the number of training data, and，("%&)
𝑤$
is the j-th weight.
("%&)

𝐾 '𝒙("%&) , 𝒙𝒋

* represents the basis function at the current time and is given by("%&) +

("%&)

𝐾 '𝒙

("%&)
, 𝒙𝒋
* = exp 5−'𝒙("%&)) − 𝒙$

𝑟$,

5

("%&)

* '𝒙("%&)) − 𝒙$

*
8

(2)

JSST2021("%&)
In the above equation, 𝑥$
and 𝑟$

are the center and radius of the j-th basis function. Learning in RBF network is to find the weight vector 𝑤 ("%&) in equation (1), which is obtained by
the following equation.
%*𝑤 ("%&) = ;𝑲("%&)+ 𝑲("%&) + 𝚲? 𝐾 ("%&)+ 𝒚(")
𝑲("%&) , 𝜦 in the above equation are expressed by the following equations, respectively.
("%&)

("%&)

("%&)

("%&)

("%&)

(3 )

("%&)

⎡𝐾 '𝒙* , 𝒙* * 𝐾 '𝒙* , 𝒙, * ⋯ 𝐾 '𝒙* , 𝒙( *⎤⎢("%&) ("%&)("%&) ("%&)("%&) ("%&) ⎥𝑲("%&) = ⎢𝐾 '𝒙, , 𝒙* * 𝐾 '𝒙, , 𝒙, * ⋯ 𝐾 '𝒙, , 𝒙( *⎥
⎢
⎥
⋮
⋮
⋱
⋮
⎢
("%&) ("%&)("%&) ("%&)("%&) ("%&) ⎥⎣𝐾 '𝒙( , 𝒙* * 𝐾 '𝒙( , 𝒙, * ⋯ 𝐾 '𝒙( , 𝒙( *⎦

(4)

(×(

𝜆 0 ⋯ 00 𝜆 ⋯ 0(5)
𝜦=MP
⋮ ⋮ ⋱ ⋮
0 0 ⋯ 𝜆 (×(In equation (5), 𝜆 represents the parameter for the weight. (𝑒. 𝑔. 𝜆 = 1.0 × 10%, )Therefore, the predicted value 𝑦" (".&) at time 𝑡 + 𝜏 is calculated as follows using the input
data at 𝒙(") the current time.
("%&)

𝑦" (".&) = ∑($)* 𝑤$("%&)

𝐾 '𝒙(") , 𝒙$*

(6)

3. RBFnetwork incorporationg affine transformation
3.1 Affine-RBF network(a) Absolute coordinate system(b)Rotated coordinate system

Fig.1Coordinate transformation for output vector using affine transformation
We focus on the direction of time series. The time series are fixed in an absolute coordinate
system, but they are converted to a rotating coordinate system using affine transformation (Fig.
1). In Fig. 1 the black dots indicate the output vector and the white dots indicate the predicted
points.
3.2 Conversion using Affine Transformation(a) Output vector conversion method(")
The output vector in the rotating coordinate system 𝑌$ ( 𝑗 = 1,2,3) shown in Fig. 1 (a) is
the output vector in the absolute coordinate system. It can be obtained using Eq. (7).

6

JSST2021t*𝑐𝑜𝑠 𝜃Z ("∗) [ = \𝑌$− 𝑠𝑖𝑛 𝜃𝑡
−𝑡0 𝑐𝑜𝑠 𝜃 − 𝑦0 𝑠𝑖𝑛 𝜃(")
c d𝑦$ e 𝑗 = 1,2, ⋯ , 𝑚
𝑡0 𝑠𝑖𝑛 𝜃 − 𝑦0 𝑐𝑜𝑠 𝜃
1

𝑠𝑖𝑛 𝜃
𝑐𝑜𝑠 𝜃

(7)

(b) Input vector conversion method(a) Absolute coordinate system(b)Rotated coordinate system

Fig.2 Coordinate transformation for input vector using affine transformation
The coordinate transformation of the K-th input vector x(jt,)K ( K = 1,2,!, k : j = 1,2,!, m) is performed as shown in Fig. 2.(𝑇1 , 𝑥1 , 𝐾) in Fig. 2 (a) shows the time in the absolute coordinate
system, and this point corresponds to the origin O in the rotating coordinate system in Fig. 2(b).
The K-th element of the j-th input vector in the rotating coordinate system is given by Eq.
(8).
Z(𝑡 − 𝑖𝜏) ∗𝑐𝑜𝑠 𝜃3("%4&)∗ [ = \− 𝑠𝑖𝑛 𝜃3𝑋$,3

𝑠𝑖𝑛 𝜃3𝑐𝑜𝑠 𝜃3𝑡 − 𝑖𝜏
−𝑡0 𝑐𝑜𝑠 𝜃3 − 𝑥0,3 𝑠𝑖𝑛 𝜃3("%4&)
c d𝑥$,3 e 𝑖 = 1,2, ⋯ , 𝑝(8)𝑡0 𝑠𝑖𝑛 𝜃3 − 𝑥0,3 𝑐𝑜𝑠 𝜃3
1

Then, short-term prediction using the RBF network is performed in the rotating coordinate
system.
When converting the predicted value 𝑌l in the rotating coordinate system to the absolutecoordinate system again, the calculation is performed using the equation (9) without using theinverse matrix of the equation (7).
𝑡∗𝑡𝑐𝑜𝑠 𝜃 − 𝑠𝑖𝑛 𝜃 𝑡0(9)m𝑦n = \cZ 𝑌 [𝑠𝑖𝑛 𝜃 𝑐𝑜𝑠 𝜃 𝑦0
1

Fig. 3 Relation between rotated and absolute coordinate system
Here, the relationship between the absolute coordinate system and the rotating coordinate
system shown in Fig. 3 is given by Eq. (10).

7

JSST2021
𝑡+𝜏
𝑐𝑜𝑠 𝜃p (".&) q = \
𝑦"
𝑠𝑖𝑛 𝜃

(𝑡 + 𝜏) ∗𝑡0
c r l (".&)∗ s𝑦0 𝑌
1− 𝑠𝑖𝑛 𝜃
𝑐𝑜𝑠 𝜃

(10)

3.3 Example of short-term prediction by RBF network using affine transformation
An example of short-term prediction by RBF network to Delayed logistic equation (Eq.11)
shown below and short-term prediction by RBF network using affine transformation is shown(Fig. 4). Here, 𝑎 = 1.3, 𝑏 = 1, 𝜏 = 1 initial value 𝜙 = 0.1 is used(4). The solid line shows
the time series data, and the white circle shows the prediction. When the white circle is on thesolid line, it means that the prediction is highly accurate. It is found from the numerical resultthat the prediction accuracy is improved by the proposed method.
𝑑𝑥(𝑡)= (𝑎 − 𝑏𝑥(𝑡 − 𝜏))𝑥(𝑡)
𝑑𝑡
2.5

2.5

RSME = 5.83 ´ 10-2

RSME = 0.78 ´ 10-2

2

2
1.5

1.5

1

x (t ) 1

x (t )

(11)

0.5

0.5

Time[s]

0
0

5

10

15

Time[s]

0
0

20

5

10

15

3

2

3

2

2

1.5

2

1.5

1

1

1

1

0

0.5

0

0.5

2

3

4

5

10

11

12

13

2

(a)Proposed approach

3

4

5

10

20

11

12

13

(b)Conventional approach

Fig.4 Illustrative example of nonlinear prediction of using Affine-RBF network

4. Coclusion
In this paper, we propose a method for short-term prediction using RBF network with affine
transformation. The conventional training data is fixed in the absolute coordinate system, andthe prediction accuracy is insufficient. By applying the proposed method to the Delayed logistic
equation, the prediction accuracy by the proposed method is improved than that of the conventional method.

References[1] He,X. and Lapedes,A.,Nonlinear modeling and prediction by successive approximation
using radial basis functions,Physica D Vol.70(1993),pp,289-301[2] Gan, M., Peng, H., and Dong, X. P., A hybrid algorithm to optimize RBF network architecture and parameters for nonlinear time series prediction, Applied Mathematical Modelling, Vol. 36, Issue 7 (2012), pp.2911-2919.
[3] Du, H. and Zhang, N., Time series prediction using evolving radial basis function networks with new encoding scheme, Neurocomputing, Vol. 71 (2008), pp.1388-1400.
[4] Aihara,Kazuyuki, Basics and Applications of Chaotic Time Series Analysis, Sangyou
Tosho(2000).(Japanese)

8JSST2021

機械学習を使用した和太鼓演奏者の判別
Discrimination of Japanese drum player using
machine learning
豊田 佳奈子 1,田村 裕 1*
Kanako TOYODA1,Hirosh TAMURA1*
1

1

中央大学大学院 理工学研究科

Graduate school of Science and Engineering,Chuo University*tamura@elect.chuo-u.ac.jp

Abstract. The AI boom is occurring due to improvements in computer functions. Machinelearning can play music automatically. However most of those are Western music. There is aprevious study of automatic Japanese drums performance. In this research, We distinguish between performance of Japanese drum experts and automatic performance .
Keywords: Computer simulation, Machine learning,Japanese drum

1. まえがき
近年は計算機の能力向上や深層学習の技術が発達し AI ブームが起きている。深層学
習では膨大なデータから AI が自ら学習し特徴を取得できるようになり、文字の生成や音楽の自動演奏を行うことが可能になった。この様な音楽演奏や楽譜の生成の歴史は長く、初めは音楽を作る為の作曲方法の知識を実装した[1]。次に音楽のデータベースを使用し[2]、その後統計的な確率モデルを基に自動的に制作[3]が可能になった。
しかしこの様にして制作される音楽のほとんどが西洋の音楽が多く日本古来の音楽の自動演奏というものはまだあまり事例がない。
埼玉県秩父地方に分布し古くから伝承されている秩父屋台囃子[4]というお囃子がある。これは秩父夜祭で演奏され、楽器の編成は大太鼓、小太鼓、笛、鉦である。秩
父屋台囃子は大太鼓と小太鼓が一緒に練習を行うものであるが、地方の過疎化問題な
どで継ぎ手が減少してしまい大太鼓の演奏者と小太鼓の演奏者が同時に練習する機会が減ってきている。ここで、大太鼓や小太鼓の演奏方法を機械学習[5]により学習させ自動演奏を可能にすることで、演奏者が個人で練習する機会を増やし伝統の保存に役立てようとする研究が存在する。先行研究[6][7]では和太鼓の自動演奏に焦点を当てていたが、出来上がった音声データが実際のお囃子に近いかを判定する際に秩父屋台囃子を知っている知っていないに関係なく複数人にアンケートを行うことでし9

JSST2021
か実現出来ていなかった。そこで本研究では、自動演奏が実際の演奏に近いのかを人間の感覚ではなく客観的に評価し判別する一環として、機械学習による小太鼓の自動演奏と熟練者の演奏した小太鼓の演奏を判別するプログラムの作成を目標とした。

2. 研究手順本研究で使用した特徴量は 1 分間での拍数を表す BPM、フーリエテンポグラムと呼ばれるノベルティファンクションに対し、短時間フーリエ変換を適用し得られた振幅スペクトルの２つである。使用した学習方法は、SVM[8]、決定木[9]、3 最近傍法[10]、ベイズ分類器[11]である。
1)SVM は Support Vector Machine の略で、入力されたデータを 2 グループに分ける為に境界線を引くことを考え、境界線に最も近いグループの異なる 2 点と境界線の距離がもっとも大きくなるように境界線を引くことを考えるアルゴリズムである。
2)決定木は、木構造を用いて入力されたデータに条件を与え 2 分割することを繰り返し、子ノードに割り当て、最終的に子ノードに存在するデータが 1 種類になるまで条件を与え続けるアルゴリズムである。
3)3 最近傍法は、新しいデータを入力しグループを割り当てる際に、新しく入力されたデータから最も近く、且つ既にグループを割り当てられたデータを 3 つ選び、選ばれたデータの中で数の多い方のグループに新しいデータを割り当てるというアルゴリズムである。
4)ベイズ分類器は、ベイズの定理を用い、データの特徴量は独立で互いに相関がないと仮定したとき、まず初めにデータが与えられたときの全ての推定の確率を求め、最も確率の高い推定結果を出力するアルゴリズムである。
BPM を特徴量としたとき、熟練者の演奏した小太鼓の音源データと、機械学習により作成した小太鼓の自動演奏の音源データを共に Python 用音楽解析モジュールである LibROSA [12]を利用し 120 秒分取り込んだ。その後２秒ごとに音源を区切りデー
タセットを作成した。フーリエテンポグラム特徴量としたとき、熟練者の演奏した小太鼓の音源データと、機械学習により作成した小太鼓の自動演奏の音源データを共に
LibROSA を利用し 120 秒分取り込んだ。その後、5 秒、10 秒と音源を区切りデータセットを作成した。

3. 研究結果及び考察
BPM を特徴量として使用した時の判別精度の結果を表 1 に記した。フーリエテンポ
グラムを特徴量として使用した時の判別精度の結果を表 2 に記した。
表 1 BPM を特徴量としたときの判別結果

特徴量

学習方法

BPM

SVM

47%

BPM

決定木

57%

BPM

3最近傍

53%

10

判別精度

JSST2021
表 2 フーリエテンポグラムを特徴量としデータセット数を変化させたときの判別結果

特徴量

学習方法

判別精度(5秒ごと区切った) 判別精度(10秒ごと区切った)フーリエテンポグラム

svm

100%

100%

フーリエテンポグラム

決定木

93%

100%

フーリエテンポグラム

3最近傍

100%

100%

フーリエテンポグラム

ベイズ分類器98%

99%

BPM を特徴量として使用したとき表 1 を見ると判別出来ている確率は 50％前後と全く判別することが出来なかった。これは使用した熟練者の音源と機械学習による自動演奏の音源では図 1 及び図 2 を見ればわかるようにテンポに違いがほとんどなく 1
分間での拍数を表す BPM では判別することが出来なかったことが原因だと考えられる。図 1 及び図 2 は横軸が時間(seconds)であり縦軸は振幅である。テンポの違いから演奏者を判別したかったので新しくフーリエテンポグラムを特徴量として使用した。このとき表２よりどの学習方法でも熟練者の演奏と機械学習による自動演奏を判別出来ている事が分かる。フーリエテンポグラムは発音のタイミングが分かりやすいように加工した信号であるノベルティファンクションに対し、短時間フーリエ変換を適用し得られた振幅スペクトルである。本研究では音源を 5 秒又は 10 秒ごとに配列に格納しそのデータにフーリエテンポグラムを適用させた。これは区切られた音声データごとに窓関数をかけ周波数分析を行うことになるため、発音のタイミングと小太
鼓演奏の音の高低を考慮することができ、特徴量を上手く作れたのではないかと考えられる。更に図 3 及び図 4 を見ると、熟練者のスペクトログラムでは 0[dB]のときの周波数は 384[Hz]であり、自動演奏のスペクトログラムは 0[dB]のときの周波数は128[dB]となっている。このことからも演奏の特徴の違いが表れている。図 3、図 4 共に横軸は時間(minutes)を表し縦軸は周波数(Hz)を表しており、色で 1 分間での拍数を表している。
表 2 を見ると曲を区切る長さを 5 秒から 10 秒にすると精度が高くなっていることが分かる。これを確認するために今までの 3 つの学習方法とは異なるベイズ分類器という学習方法を使用し、この学習方法でも判別精度を上げることが出来た。以上から、入力した音源を短時間フーリエ変換する際に 5 秒ずつ区切ったものよりも 10 秒ずつ
区切ったものの方が音声信号の情報量が大きくなり特徴量として良いものになったのではないかと考えられえる。

図 1 時間ごとの熟練者の演奏の振幅を表す波形

図 2 時間ごとの機械学習による自動演奏の振幅を表す波形

11

JSST2021

図 3 熟練者の演奏のスペクトログラム

図 4 自動演奏のスペクトログラム

4. 今後の課題
今後の課題としては以下の２点が考えられる。1 つ目は、今回の研究では秩父屋台
囃子の小太鼓の音声データが熟練者によって演奏されたものか、機械学習による自動演奏かを判別することしかできない。よって機械学習による演奏が上達しているかを判別することは実現出来ていないので今後の課題である。2 つ目は、今回は秩父屋台
囃子の小太鼓の演奏に限定して研究を行ったが他の打楽器でも違いを判別できるか検証することも今後の課題である。

参考文献[1] Sundberg, J. , Askenfelt, A. and Fryden, L.: Musical Performance, A Synthesis-by-RuleApproach, Computer
Music Journal, Vol.7, pp.37-43(1983).
[2] Cope, D.: Machine Models of Music, MIR Press, pp.403-425(1992).
[3] Pachet, F. and Roy, P.: Markov Constraints: Steerable Generation of Markov Sequences,Constraints, Vol.16,
No.2, pp.148-172(2011).
[4] www.asahi-net.or.jp/~qj4t-nkmr/yataibayashi2/yataibayashi.htm[5] https://www.codexa.net/what-is-machine-learning/[6] 佐藤信太郎 "和楽器における機械学習を用いた演奏支援システムについて" 中央大学理工学研究科電気電子情報
通信工学専攻修士論文 2019[7] 岡本和巳 "機械学習を用いた和太鼓の自動演奏について"中央大学理工学部電気電子情報通信工学科
2019[8] https://avinton.com/academy/svm[9] https://qiita.com/3000manjPY/items/ef7495960f472ec14377[10] https://qiita.com/yshi12/items/26771139672d40a0be32#:[11] https://qiita.com/aflc/items/13fe52243c35d3b678b0[12] https://www.wizard-notes.com/entry/music-analysis/librosa

12

卒業論文

JSST 2021

移動音の移動方向知覚に関する研究-適応法による検討Perception of Direction of Moving Sounds.
-Study on Adaptive Method瀧澤 哲 1,* , 工藤 彰洋 1
Tetsu Takizawa1,* , Akihiro Kudo1
1 苫小牧工業高等専門学校
1 National Institute of Technology, Tomakomai College∗ tt20805@tomakomai.kosen-ac.jp

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

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

13

JSST 2021

第 5 世代移動通信システムにおけるエッジコンピューティング用サーバの設置方法について
A Method for setting up Servers for Edge Computing
in Fifth-Generation Mobile Communication Systems
江原 亮太 1, 田村 裕 1*
Ryota EHARA1, Hiroshi TAMURA1*

1
1

中央大学大学院 理工学研究科

Gladuate School of Science and Engineering, Chuo University*

tamura@elect.chuo-u.ac.jp

Abstract. While the 5th generation mobile communication system (5G) has improved communication speed and capacity compared to the conventional communication systems, it also has
the problem of reduced communication range and obstacle tolerance. To solve this problem, amethod called mobile edge computing (MEC) has been proposed, in which a Simplified serveris installed close to the terminal to perform some of the processing on its behalf, thereby reducing the burden on the main server and shortening the communication distance at the sametime. In this paper, we will examine the location of a simple server that enables efficient communication under the assumption that the server is installed in a base station.
Keywords: Fifth Generation Mobile Communication Systems, Mobile Edge Computing

1. まえがき
従来の移動通信システム(4G や LTE と呼ばれるシステム)より通信能力を向上させた
第 5 世代移動通信システム(Fifth Generation Mobile Communication Systems:以降 5Gと表記する)が現在普及に向けて整備が進められている。5G は通信速度や通信容量の向上を達成した反面、通信範囲の縮小や直進性の強化による対障害物耐性の低下、通信量の肥大化によるメインサーバの処理能力の限界による処理作業の遅延といった問題が生じている[1]。前者 2 点の問題については中継用の基地局を増やすことで対処している。一方後者の問題を解決するにあたり、端末－メインサーバ間の通信路上に簡易的な処理用サーバを設置して一部の処理を肩代わりさせるモバイルエッジコンピューティング(Mobile Edge Computing System:MEC)という技術の導入が考えられている[2]。

17JSST 2021
本稿では、この簡易サーバ(以降 MEC サーバと表記する)の配置について、中継用の基地局の一部に MEC サーバを設置する条件下で、通信効率の良いサーバの設置場所を検討していく。

2. シミュレーション
高速通信を保持するためには、基地局から一定の距離以内に MEC サーバを設置する必要がある。この時必要なサーバの数が少なくなるように設置方法を検討する。
1×1[km2]の 2 次元空間をシミュレーションエリアとし、エリア内に母点をランダムに配置する(図 1)。この母点を基地局の設置箇所とし、エリア内のすべての基地局
が MEC サーバを仲介して相互に高速通信を行える状態が成立するよう、MEC サー
バを基地局上に設置していく。5G での高速通信を成立させる条件として、静止時に通信距離が約 0.3[km]の状態で成立している[3]ことから、本稿では半径 0.25[km]圏内
で高速通信が成立するという条件のもと MEC サーバを設置する。また、MEC の並列処理の上限を直接接続する基地局 5 基分と仮定してシミュレーションを行う。
シミュレーションの初期状態として、MEC サーバが設置されておらず、最初のサーバ設置場所を任意の基地局から 1 基選ぶ場合(パターン①)と、あらかじめなるべく均等な配置になるよういくつかの MEC サーバを設置した状態を考え、後者について4 等分した場合(パターン②)と 9 等分した場合(パターン③)の 2 つのパターンで初期の MEC サーバの配置を考える。
まず、設置された MEC サーバの高速通信圏を確認し、サーバと接続されていない基地局から優先的にサーバの接続上限まで接続していく。その後エリア内の全てのMEC サーバの高速通信圏内にある基地局 (図 2 の黒円内の母点)からサーバ未接続の基地局(マークのついてない母点)を優先してランダムに 1 基選択し、新たな MEC サーバを設置する。この作業を繰り返し、エリア内のどの MEC サーバにも接続していない基地局が存在しないかつすべての MEC サーバが高速通信のみで相互接続できる状態になった時にシミュレーションを終了する。
これらのシミュレーションを 10 通りの基地局の配置で実行し、MEC サーバの設置数の平均をとる。また、基地局数を変化させた場合について同様にシミュレーションを行う。
シミュレーションの設定を表 1 にまとめる。

図 1．シミュレーションエリアの例

図 2．MEC サーバ(赤円)への接続例(基地局 30 基)(接続済み基地局:青円)

18

JSST 2021
表 1．シミュレーション設定
シミュレーションエリア

1[km]×1[km]の2次元空間

基地局(母点)配置

ランダム

高速通信圏

半径0.25[km]圏内

MEC－基地局間接続条件

MECの高速通信可能領域内

MECサーバの接続上限

基地局5基分(サーバを設置した基地局を含む)

MECの初期配置パターン①0基、②5基(4等分＋エリア中心)、③9基(9等分)MEC設置割合算出方法(MEC設置数/基地局数)×100[%](下限20%)

3. 結果と考察
表 1 の設定を適用し、基地局数が 30 基の場合および 50 基の場合についてシミュレーションを行った。シミュレーション結果の一例を図 3～図 6 に示す。また、シミュ
レーション結果について、基地局の総数に対する MEC サーバの設置割合の平均を表2 に示す。
基地局数が 30 基の場合、あらかじめ MEC を設置するパターン②や③より MEC サーバを設置していない状態から開始したパターン①のほうが MEC サーバの総設置割合が明らかに少なくなる傾向がみられる。基地局数が少ないとき、配置状況次第では基地局が密集している空間や反対に基地局がほとんどない空間が生じる。均等配置では基地局の配置状況に関係なくエリア内になるべく等分に MEC サーバを設置した状
態を初期状態とするため、密集地帯に設置されたサーバがカバーしきれずに追加設置量が多くなったり、基地局間距離が離れているために MEC サーバ同士の高速通信を成立させるためにサーバを多数設置する必要ができてしまうことが要因となり、結果として MEC サーバの総設置割合が多くなる傾向がみられるものと考えられる。
基地局数が 50 基の場合、傾向は 30 基の場合に類似しているが、MEC サーバの総設置割合の差が小さくなっている。これは基地局数が増えたことにより、パターン①の場合でも必要な MEC サーバの設置数が増え、あらかじめ設置しているパターン②
や③の状況に近づいたためと考えられる。

図 3．パターン①(30 基)の例

図 4．パターン②(50 基)の例
(赤円は初期配置上の MEC サーバ)

19

JSST 2021

図 6．パターン③(50 基)の例

図 5．パターン③(30 基)の例

表 2．シミュレーション結果

シミュレーション
MEC総配置割合

パターン①

パターン②

パターン③

30基

50基

30基

50基

30基

50基

34.30%

23.20%

59.30%

25.60%

47.00%

33.20%

4. まとめと今後の課題
本稿では、MEC サーバを用いた高速通信を実現する際、MEC サーバの設置数を最小化するためにはパターン②や③が必ずしも有効な手法ではないという結果が得られた。
今回のシミュレーションでは、エリア内に通信の妨げとなりうるような障害物が配置されていない平面空間を舞台としたが、現実空間により近いシミュレーションを行
う上では電波を遮断しうる障害物が存在する空間でのシミュレーションが必要であると考えられる。また、基地局以外のところに MEC サーバを配置していく場合についての設置手順の考案も今後の課題として考慮していく必要がある。

参考文献
[1]

[2]

[3]

5G の商用化から普及に向けての課題：株式会社 日立ソリューションズ・クリエイト
https://www.hitachi-solutions-create.co.jp/column/mobile/5g-task.html
最終閲覧日：2021 年 7 月 14 日
5G の技術を整理する 当初はなぜ高速大容量通信しか使えないのか？：5G ビジ
ネスの神髄に迫る - ITmedia Mobile
https://www.itmedia.co.jp/mobile/articles/2004/24/news083_2.html
最終閲覧日：2021 年 7 月 14 日
岸山祥久, 原田篤, 余波, ベンジャブールアナス, 奥村幸彦:” ミリ波を用いた超高速・長距離伝送の５G 屋外実験”, NTT DOCOMO テクニカルジャーナル Vol.26 No.1
PP.25-32, 2018.

20

JSST 2021

グラフ畳み込みネットワークによる数独の解法
Solving Sudoku
with Graph Convolutional Networks
石井 宏一郎 1, 田村 裕 1*
Koichiro Ishii1, Hiroshi Tamura1*
1

中央大学大学院 理工学研究科*tamura@elect.chuo-u.ac.jp

Abstract. Sudoku is a puzzle type in which a 9 × 9 grid is filled in with numbers, based on numbers that have already been filled in. Graph convolutional networks (GCN) have recently become
a popular method for deep learning using graph representations. This paper considers Sudoku as agraph structure and tries to these puzzles using GCN. In solving Sudoku puzzles, individuals use
the candidate numbers in the blank fields as clues. With a GCN, these candidates can be regarded
as features for learning. Therefore, we thought this structure would more easily reveal Sudoku’s
features for learning and how they operate, making the learning process more closely resemble theway humans solve puzzles. Additionally, Sudoku puzzles have recently been created with different rules, and this GCN model is sufficiently flexible to be applied to these new puzzles. We prepared some learning models with different structures and feature values. In the paper’s final section, we compared and discussed differences in these models’ learning accuracy.
Keywords: Machine learning, Sudoku, GCN

1. まえがきディープラーニングによる機械学習は技術の発展と共に大きく進化を遂げたと同時に、世間一般にも広く知られるようになった。この要因の一つとして、画像処理分
野への機械学習の適用が大きな成功を収めたことが挙げられる。この分野では順伝搬型ニューラルネットワークの一つ、畳み込みニューラルネットワーク(CNN; Convolutional Neural Network)が主流に使用されている。ここで用いられる畳み込み演算は、
入力データに対して同じ次元を持ったフィルタをずらしながら積和計算するため、データの形状を維持したまま学習が可能なことが特徴である。しかし、この演算は格子状などの一定の規則性が入力データの形状に必要とするために、任意の形をしたデータに対しては適用できなかった。そこで任意の構造を持ったデータ、主にグラフ構造に対しても畳み込みニューラルネットワークと同様の処理によって学習する方法がこれまでに様々考案され、その一つがグラフ畳み込みネットワーク(GCN; Graph

21

JSST 2021
Convolutional Network)である[1]。
本論文ではグラフ畳み込みネットワークを利用し、数独[2][3]の解答を試みる。数独の解答への畳み込みニューラルネットワーク(CNN)の適用はこれまでいくつか存在し[4][5]、99%の精度を出している。この方法ではフィルタを用いて数独の問題の各行、各列、各ブロックでの特徴を取り出すことで学習する。この時フィルタは複数個用意され、この数によって精度が変動するため、各区画の特徴を正しく把握することが難しい。そこでグラフ畳み込みネットワークを用いることによって各マスの接続関係から得られる各マスの特徴を数独に適用することで、人間が解く際に参考にする「候補」と特徴量を対応づけることができ、特徴を理解しやすく、また解き方を人間らしさに近づけることができると考えた。

2. グラフ畳み込みネットワーク
ノードとノード間を結ぶエッジから構成される構造をグラフ構造と呼ぶ。グラフはノードごとの接続関係に規則性がないため、通常の畳み込み演算をそのまま適用することができない。そこで接続しているノードからと自身のみの情報を足し合わせた後、重みを乗じることで自身の情報を更新するという操作により、畳み込み演算を定義している。ノード i における畳み込み演算は以下の式で表される[1]。
(𝑙+1)

ℎ𝑖

1
(𝑙) (𝑙)(𝑙) (𝑙)𝑊1 ℎ𝑗 + 𝑊0 ℎ𝑖 )𝑐𝑘 𝑖,𝑗= 𝜎( ∑𝑗∈𝑁𝑖

ここで𝑐𝑖,𝑗 は正規化をするための定数、𝑁𝑖𝑘 はノード i から k 以下の距離でつながっている近傍ノードの集合、𝑊 (𝑙) はネットワークの l 層における学習重み行列、σ(∙)はネットワークの出力を決定する活性化関数を表す。
また、この畳み込み処理の概略図を以下図 1 に示す。各点において隣接している点から矢印の方向に情報を取り込む。最後に活性化関数に通して１層の処理が終了する。
本論文では活性化関数として、０以下なら「０」、０よりも上であれば「入力値と同じ値」を返す関数 ReLU を用いている。

図１

グラフ畳み込み処理

22

JSST 2021

3. 数独の解法への適用
人間が実際に数独を解く際、それぞれのマスに入る数字の候補を絞っていく。この解き方と同様に、初期配置から各マスの候補をベクトルに対応させたものを各点の特徴とし、グラフ畳み込みネットワークを構成することで、数独の学習を試みた。すでに最初に数字が配置されているマスは情報を更新する必要はないが、そのままの構造であるとすべての点（マス）が情報を更新する対象となるため、問題の空欄のみが更新し続けるような構造にする。ネットワークの出力と正解ラベルとの差を交差エント
ロピー誤差として表すと、繰り返し回数(横軸)による誤差(縦軸)の推移は以下のようになった。また、ネットワーク層は３層である。

図２ ネットワークの学習

以上の結果よりエントロピー誤差は全体的には減少しているが振動は続いてしまっており、問題によって精度に大きなブレができてしまっていると考えられる。
次に、実際にはどの程度の正解率になるのかを調べるために、教師データとテストデータからそれぞれ１問を選び、その問題の正解率が繰り返し回数によってどのように変化するかを調べた。縦軸は正解率、横軸はエポック数をそれぞれ表す。

図３

同じ問題の１エポックごとの正解率

23

JSST 2021
図より、およそ８割の正解率を保ってはいるが、すべてのマスを正解させることはできていない。人間が解くときは１マス１マスを確実に埋めていく、この方法を真似して、ネットワーク出力の確率が高いマスから１つずつ埋めていくようなループ構造にしてみたが、結果は図３と同じであった。結果としてこの構造では完全に一問を解くことは難しい。

4. まとめと今後の課題
数独をそのままグラフ畳み込みネットワークに適用させたときの学習は、結果としてまだ精度の面に難が残る結果となった。数独はルールとして数字に意味を持たず、各数
字の場所の入れ替えが可能なほか、記号などで代替することができる。今回はこのような性質を学習が不十分であったと考えられる。学習データの難易度を非常に簡単にする
ことや、候補の表し方をベクトルでなく、１つの数にまとめることで、点や辺の特徴をさらに特徴づけたものとすることなどをすることにより、この性質を十分に学習させる
改善が今後の課題となる。
数独の問題の難易度は、問題の空欄の数によってある程度分類されているが、この空欄の数と難易度は一般的には比例しないため、製作者によって決定されることが多い。
このように数独の難易度はいまだ明確に定義できていない。ニューラルネットワークにより難易度を計算する手法なども提案されているが[6]、特徴の人間らしさをもつ GCN に
よって数独の解法、テクニックを何らかのフィルタとして表現することによって、将来解法の難易度などの人間が理解できる根拠による難易度の定義ができるのではと考えて
いる。

参考文献
[1]立花誠人 村田剛志: 構造特徴とグラフ畳み込みを用いたネットワークの半教師あり学習 人工知能学会論文誌 34(5), B-IC2_1-8, 2019[2]Tom Davis: “The Mathematics of Sudoku | Mathematical Circles Topics”
http://www.geometer.org/mathcircles/ (2021 年閲覧)
[3]数独 - Wikipedia
https://ja.wikipedia.org/wiki/%E6%95%B0%E7%8B%AC (2021 年閲覧)[4] Shiva Verma: “Solving Sudoku with Convolution Neural Network | Keras” (2019)https://towardsdatascience.com/solving-sudoku-with-convolution-neural-network-ker
as-655ba4be3b11 (2021 年閲覧)[5] GitHub, “nnethercott /Sudoku-CNN”
https://github.com/nnethercott/Sudoku-CNN (2021 年閲覧)[6]天野秀亮 篠埜功 杉本徹: ニューラルネットワークによる数独の難易度判定手法の提案 信学技報, vol. 112, no. 319, AI2012-17, pp. 13-18, 2012 年 11 月.

24

JSST 2021

無人航空機を利用した中山間地域における配達についてOn delivery methods using UAV in hilly and mountainous areas
増子 裕斗 1*, 田村 裕 2
Yuto MASHIKO1*, Hiroshi TAMURA2

1
1

中央大学大学院 理工学研究科

Graduate School of Science and Engineering, Chuo University*tamura@elect.chuo-u.ac.jp

Abstract. UAV(drones) are growing rapidly and expectations for drones are rising for a varietyof reasons. This paper shows the usefulness of delivery using drones. The result of the simulation shows that it was found that four drones are required to work for one car. In addition, thefastest delivery can be complered by flying the drone in parallel with the car and delivering it.
Keywords: Computer simulation, Drone, Delivery Service

1. まえがき
近年、特に中山間地域において、「買い物難民」と呼ばれる自身で身体的理由や経済的理由等により買い物に行くことができず、日用品を買うことが困難となる人々及びその現象が問題となっている。この問題の解決策の 1 つとして流通業の強化がある。その中でも商品そのものが家の玄関までやってくる「宅配サービス」について本稿では考える。宅配業の者が利用者が頼んだ商品を届けてくれるので、利用者本人は家から出ることなく食料品や日用品などの生活必需品のみでなく、娯楽品なども手に入れることが出来る。
しかし、昨今のインターネットの普及につれインターネットショッピング等の利用者が増え、宅配業車にかかる負担が以前よりはるかに増している。その結果、宅配業
者は 1 日に 150~200 個の配達を行うことになっている。
そこで、本稿では人による配達ではなく、小型無人航空機（ドローン）を用いた荷25

JSST 2021
物配達について考える。現在ドローン市場は急速に成長しており、少子高齢化による慢性的な労働力不足や人件費高騰による省力化・無人化の促進と低コスト化などの理由によりドローンへの期待が高まっている[1]。実際に、2021 年 1 月に楽天が三重県にある離島への配送サービスを行っている[2]。本項では、まずドローン 1 台による配達に必要な時間と自動車 1 台による配達に必要な時間とを比較して、ドローン 1 台分の働きについて考察する。さらに、ドローン単独だけでなく、自動車とドローンとを併用する形で配達を行うときについても考え、その場合に必要な時間を比較すること
で、ドローンの有用性ついて検討を行う。また、それぞれの場合における環境への影
響についても考察する。

2. シミュレーション
今回のシミュレーションにおけるモデルは 10[km]×10[km]の正方形エリアとする。
その中心の座標は(5000,5000)であり、エリアの四隅の差表を(0,0)、(10000,0)、(0,10000)、(10000,10000)とする。
このエリア内において、中心に配送拠点を置き、配達利用者の家を 100 軒ランダムに配置する。このエリアについて概形を図 2-1 に示す。利用者は場合によりドローンによる配達なのか自動車による配達なのかを決定する。場合については、ドローン 1台による配達を行う場合、自動車 1 台における配達を行う場合、自動車とドローンを
各 1 台ずつ用いた場合の大きく分けて 3 つの場合について考える。自動車とドローンとを用いた場合のそれぞれの動きについて図 2-1 に追記し示す。
ドローンの速度については、実際に岩手県において実験に使用されたドローンをもとに 55[km/h]として考える[3]。また、ドローンの配達については、基本的に目的の配達が終了したら中央の配送拠点に戻り次の配達分の荷物を積むこととする。ただし、自動車とドローンを併用した場合については、ドローンは中央の配送拠点に戻らず、同エリア内を走行している自動車のもとへ飛行し、そこで次の配達分の荷物を受け取るものとする。また、自動車については、基本的に時速 60[km]で走行するものとし、ランダムで一時停止や信号待ち等を考慮した待機時間を考慮してシミュレーションを行うこととする。また、荷物の積み下ろしや受け渡しの時間は考慮しないものとする。
（0,10000）(10000,10000)(0,10000）(10000,10000)

10
km

(0,0)
:配送拠点

10
km

(5000,5000)(5000,5000)

10km(0,10000)
(0,0)
10km
(0,10000):サービス利用者の家:ドローンの経路:自動車の経路
図 2-1 シミュレーションモデル

26

JSST 2021

3. 結果と考察
二次元エリア下におけるシミュレーションの結果を表 3-1 に示す。
まず、ドローン及び自動車について、各々一台のみを用いた配達について考える。
ドローン一台のみでの配達を方法１、ドローン一台のみでの配達だが、一度に 2 件分の荷物を配達できるものとする配達を方法２、自動車一台の配達を方法３とする。本稿では、ドローンは一度に 1 件または 2 件分の荷物しか積むことができず、その都度荷物を積みに配送拠点に戻らなくてはならないため、サービス利用者の家から家へと配達を行うことが出来る自動車に比べ 2~3.5 倍の時間がかかってしまうと考られる。
つまり、自動車一台分の働きをするためには 2~4 台のドローンが必要となる。
次に、ドローンと自動車を一台ずつ利用した場合の方法４について考える。100 件
全ての配達を終える時間は最も速くなっていることがわかる。特にドローンによる配達については、ドローン 1 台による配達にかかる時間と比べ 4 分の 1 倍と大幅に短縮されていることがわかる。これは、単純に自動車との併用による配達件数の減少もあるがそれと同時にドローン 1 台による配達の時とは異なり、ドローンは毎回の配達時にエリア中央の配送拠点に戻らずに、エリア内を走行する自動車の元へと飛行するため、このような結果になったと考られる。ドローンによる配達を行う場合、どうしても一度に多量の荷物を積むことが不可能であるため、往復の行程が必須になってしまう。配達にかかる時間を少なくするということは、この往復の工程に必要な時間を少なくするということであるので、今回のように自動車に並走する形で配達を行う
ことで、配達全体に必要な時間を少なくすることが出来ている。ドローンを用いた配達を行う場合、必要な時間を少なくするにはドローンが何も持たずに飛行している時間を短くすることが重要だと考えられる。また、配達時間を減少することで、ドライ
バーへの負担を減らすことができるが、ドローンを操縦する人間を別に配置する必要がある。
本条件のもとで、ドローン一台あたりの二酸化炭素の排出量は自動車一台の二酸化炭素の排出量と比べ、配達件数が同数の場合ドローンの方が約 2[kg]程少ない[4]。二酸化炭素の排出量を一日の配達で約 2[kg]減少できるとすると年間で 700[kg]以上減少することができる。
表 3-1.シミュレーション結果
配達にかかった時間
シミュレーション条件

二次元エリア
ドローン

自動車

方法１

約 14 時間―

方法２

約 9 時間―方法３―

約 4 時間方法４

約 3.5 時間

約 3.5 時間27

JSST 2021

4. まとめと今後の課題
本論文では、ドローンを利用した効率の良い配達の仕方について検討を行った。ドローンをうまく活用することで、配達全体にかかる時間を短くできることがわかった。また、ドローンを使い配達時間を減少させることで、CO2 排出量についても減少できることがわっかたので、ドローンによる配達は環境に優しいという側面も持つといえる。
今後の課題として、地形による変化であったり、家の配置に偏りがあったりと、今回シミュレーションを行った場合以外の点にも考慮しなくてはならない。それを考慮したシミュレーションを行うことでより精度の高い結果を得ることができると考えている。

参考文献[1] 野波健蔵 「「空の産業革命」をもたらすドローンの課題と展望」
計測と制御 56 巻 1 号 P2 2017 年
https://www.jstage.jst.go.jp/article/sicejl/56/1/56_2/_article/-char/ja/[2] 「楽天、自動飛行ドローンによる配送サービスを離島で提供」楽天株式会社
2021 年 1 月 6 日
https://corp.rakuten.co.jp/news/update/2021/0106_01.html[3] 「DJI 社製「QS8」を使用した、孤立病院への救援物資輸送に成功！」株式会社
シーカー 2018 年 12 月 10 日
https://prtimes.jp/main/html/rd/p/000000019.000025357.html[4] Jennifer Langston「Drone vs. truck deliveries: Which create less carbon pollution?」
UNIVERSITY of WASHINGTON 2017 年 3 月 30 日
Drone vs. truck deliveries: Which create less carbon pollution? | UW News (washington.edu)

28

JSST 2021

複合現実技術を用いた放射線可視化システムの開発
Development of visualization system for radiation
by mixed reality technology
飯田 聖人 1, 工藤 颯太郎 1, 河原 梨花 2, 山口 克彦 3, 齋藤 誠紀 1*
Masato Iida1, Sotaro Kudo1, Rika Kawahara2, Katsuhiko Yamaguchi3, Seiki Saito1*

1

山形大学 工学部 情報・エレクトロニクス学科
福島大学院 共生システム理工学研究科 環境放射能学専攻
3
福島大学 理工学群 共生システム理工学類
1
Department of Informatics and Electronics, Faculty of Engineering, Yamagata University
2
Major in Environmental Radioactivity, Faculty of Symbiotic Systems Science, Fukushima University
3
Faculty of Symbiotic Systems Science, Cluster of Science and Technology, Fukushima University*saitos@yz.yamagata-u.ac.jp
2

Abstract. The accident at the Fukushima Daiichi Nuclear Power Plant on March 11, 2011, releaseda large amount of radioactive materials. It is difficult to work safely only by measuring radiation
doses. Moreover, the visualization of the radiation enhances the safety of the work. Therefore, wedevelop a system that display radiation trajectories and doses by mixed reality technology. In thissystem, HoloLens2 scans QR codes which contain the information of radiation source. The trajectories of radiation are calculated by a Monte Carlo code PHITS to display on the HoloLens2 screen.
Keywords: Radiation, Mixed reality, HoloLens

1. 研究背景・目的2011 年 3 月 11 日に発生した福島第一原子力発電所事故により，大量の放射性物質が放出された．その影響により，バリケードなど物理的な防護措置が実施され，避難が求められている区域が現在も存在する．居住可能な区域が増えるように除染作業が実施されている．
除染作業は線量計測を実施しながら行われる．しかし，放射線を直接目で見ることはで

29

JSST 2021
きないため，安全な場所を把握しながら作業することが難しい．そこで，本研究では，複合現実（Mixed Reality: MR）技術を用いて除染作業を補助するシステムを開発する．仮想空間内に放射線の軌跡や線量を表示し，放射線源や線量を目視で確認できる除染・廃炉作業補助 MR システムの構築を目指す．

2. 放射線可視化システムの概要
本研究では HoloLens2 という Microsoft 社の MR ゴーグルを用いる．開発している MR
システムは主に HoloLens2，計算サーバからなる．計算サーバは，PHITS (Particle and Heavy
Ion Transport code System) [1] と PHITS 制御コードの 2 つの要素から構成される．PHITS
は「あらゆる物質の中で放射線挙動を核反応モデルや核データなどを用いて模擬するモ
ンテカルロ計算コード」である．α 線，β 線，γ 線，中性子線などの放射線のシミュレー
ションが可能である．放射線可視化システムは以下のように動作する．
① HoloLens2 の環境認識カメラを用いて線源，遮蔽物の種類，形状，エネルギーが記録された QR コードを読み取る．
② 入手した情報を計算サーバに送る．
③ 入手した情報から計算サーバで PHITS 制御コードを実行して PHITS 実行に必要な入力ファイルを作成する．
④ 計算サーバで PHITS を実行して線源から放射された放射線の軌跡を計算する．
⑤ 計算結果をダウンロードし，HoloLens2 の画面に表示する．
本システムで可視化した中性子線の様子を図 1 に示す．円筒形の線源から放出される
中性子線が空気中を進む様子が確認できる．
現在は QR コードを用いて線源，遮蔽物の情報を読み取っているが，将来は計測した線量をもとに線源を特定 [2] することを検討している．また，周辺環境の三次元情報を計算に取り入れ遮蔽物によって放射線が遮蔽される様子も計算することを予定している．
本研究における線源，遮蔽物等の情報の送信とサーバで実行されたシミュレーション結果の受信を第 5 世代移動通信システム(5G)を用いて高速に行うことを目指す．

図 1 HoloLens2 で可視化した円筒形の線源から放射される中性子線の様子．

30

JSST 2021

参考文献
[1] T. Sato, Y. Iwamoto, S. Hashimoto, T. Ogawa, T. Furuta, S. Abe, T. Kai, P. E. Tsai, N.
Matsuda, H. Iwase, N. Shigyo, L. Sihver, K. Niita: Features of Particle and Heavy Ion
Transport code System (PHITS) version 3.02, J. Nucl. Sci. Technol. 55:5-6 (2018), 684690.
[2] T. Uemura, K. Yamaguchi: Estimation of radiation source distribution using ma-chine learning with γ ray energy spectra, J. Adv. Simulat. Sci. Eng. 7:1 (2020),
71-81.

31

JSST 2021

Computation of Micro-particle Dynamics in a
Filtration Process with the Coulomb Interaction
Nobuyuki Hirooka1,* and Satoru Ushijima2
1 Department of Civil and Earth Resources Engineering,Doctoral Program,Kyoto University
2 Academic Center for Computing and Media Studies (ACCMS), Kyoto University∗ hirooka.nobuyuki.85e@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.14.96

32

JSST2021

UAV assisted Survivor Search Scheme with Survivor
Density Distribution Prediction
Chen Haosong1*, Kazutoshi Yoshii2,4, Shigeru Shimamoto3
Department of Computer Science and Engineering, School of Fundamental Science and Engi-

1

neering, Waseda University*hchen53@ruri.waseda.jp, kazutoshi@suou.waseda.jp,

shima@waseda.jp

Abstract. UAVs, which stands for the unmanned aerial vehicle, are being used more and morein disaster search and rescue (SAR) missions due to their fast response time and wide field ofview. However, using UAVs in disasters also comes with issues of time and energy constraints.
In this paper, we addressed the problem by utilizing the call detail records(CDR) to extract thecellular activity distribution beforehand, and use it in machine learning to predict the future
cellular activity distribution, which can be used by the rescue workers to understand what thesurvivor distribution is like, and plan an optimal UAV SAR mission based on it.
Keywords: UAV, Machine Learning, Call Detail Records, Flight Simulation

1. Introduction
UAV, which stands for the unmanned aerial vehicle, has been used more and more in the civil
field. And one of such use cases is using the UAV for search and rescue missions during natural
disasters.
However, one of the main challenges of using UAVs in the search and rescue mission is theenergy limitation [1]. Currently, most UAVs use batteries as their energy source. On the other
hand, a search and rescue mission is a time-critical mission. Less search time can mean a higher
probability of finding survivors [2].
Given the above challenges, this paper proposed a new method to conduct UAV’s search andrescue mission by using the data in the call detail records to generate a survivor density distribution map before the mission start and thus giving rescue workers a clearer image of the disaster-stricken area to plan to an optimal search and rescue mission under a limited time frameand energy source.
We also did UAV flight simulation tests in our research to find out the best search pattern.

2. Call Detail Record and datasheet description
Call Detail Record(CDR) is the data record that is generated by telephone switch and contains
the important detail of that particular telecommunication activity, such as the approximate location of the activity and the destination of the activity [3]. In other words, each time you make36

JSST2021
a phone call, send a message, or browse a webpage, a cell tower will be connected. And the
database of the mobile network will record a metadata, which is the call detail record. Knowing
the approximate location of each mobile communication activity and the universal adoption ofsmartphones means the data generated in the mobile communication network may be one of
the best ways to estimate survivor distribution patterns after natural disasters.
The CDR dataset we used in this research was obtained from Telecom Italia which has been
made public [4]. This dataset contains cellular activities of the city of Milano over 2 months oftime.
We aggregate all the cellular activities based on their associated location in the SQL database,to calculate the total number of cellular activities in each grid cell. Then we used QGIS and the
OpenStreetMap generated heat map of this result. Fig.1 shows the mentioned heat map and themetric of the heat map.

3. Machine learning on the data
Although we can extract the cellular activity distribution map during the normal time, thismight not be the case for the time of natural disaster. It’s possible that telecommunication infrastructure can be damaged in disasters [5]. And thus makes it hard or impossible to extract
data from it during the time of disaster. We proposed a method to use past cellular activity
distribution as training material to train a machine learning model, which will be able to predictthe cellular activity distribution during the time of natural disaster. The rescue worker or the
UAV operator can still use it as a reference to plan an optimal UAV search and rescue mission
even if the telecommunication infrastructure is damaged.
Since we extract the cellar activity distribution from over 60 days of CDR and with dates, theFig. 1. The visualized heatmap

37

JSST2021
data we used here is the time series data, which is a series of data points with timestamps. In
order to use machine learning on time series data, we need to convert it into supervised data
first. The method we used to convert the data is called the sliding window, which used previous
time steps to predict the next time step. Sliding window can convert a sequence of interrelated
data points into supervised data.
We used 6 different machine learning models on our data to predict future cellular activity
distribution during disasters, which are baseline model, ARMA, random forest regressor,
XGBoost, LSTM and 1 dimensional CNN. We also calculate the MSE(mean squared error)
value of each model for evaluation purpose.
Fig. 2 compared the MSE value result of all models we used in our research. For the baselinemodel or the naive model, we used it as the evaluation reference for other models. The baselinemodel has the MSE value of 167, and other models’ error values are all lower than it, whichmeans it’s reasonable to use other models in our research. Since the data we used is stationary,
ARMA also had a good performance with the MSE of 112, because ARMA is good for predicting stationary data. Random Forest and XGboost also performed very good, and outperformed
the ARMA model in our case with the MSE of 67 and 52. Random Forest uses bagging in the
training process which reduces the variance [6]. We believed this can be one of the factors that
contributed to its good performance in our case. XGboost outperformed Random Forest slightly
in our case. For LSTM model, we used the typical vanilla LSTM model, which consists of asingle hidden layer and an output layer. We set the model with 100 LSTM units in the hidden
layer to do the prediction. After training the model with our data, we got the second best result
among all models in terms of MSE value with the MSE of 25. For the CNN model, since ourdata is time series data which is a sequence of input. We used the 1 dimensional CNN in which
the kernel only moves in 1 direction and is usually used for time series prediction. Our 1D CNN
consists of 1 convolutional layer, 1 maxpooling layer, 1 flatten layer, 1 fully connected layer
and an output layer. The 1D CNN has a similar performance compares to the LSTM modelwith the MSE value of 18.

4. UAV Flight Simulation
We also did simulated UAV flights to find out the optimal search pattern for UAV’s search andrescue operation during disaster. We simulated 3 search patterns, which are parallel track in
Fig.3, creeping line in Fig.4 and expanding square in Fig.5. These are common patterns usedby rescue teams. We used the UAV simulation software SITL(software in the loop) simulator

Fig. 2. Result Comparison

38

JSST2021

Fig. 3. Parallel Track

Fig. 4. Creeping line

1

Fig . 5. Expanding Square

Fig . 6. Result Comparison

to do the flight simulation. For the simulation parameters, we set the altitude to 100 meters, thesearching area to a 0.6 kilometers area, the UAV type to quadcopter, and the ground speed to 5
meters per second. We showed the result in Fig. 6. Parallel Track has the best performance withthe time consumption of 1400 seconds, and creeping line has the worst performance with thetime consumption of 2100 seconds.
1

References[1] A. L. Adams et al., “Search is a time-critical event: when search and rescue missions
may become futile,” Wilderness Environ. Med., vol. 18, no. 2, pp. 95–101, Summer 2007.
[2] J. Steenbruggen, E. Tranos, and P. Nijkamp, “Data from mobile phone operators: A tool
for smarter cities?,” Telecomm. Policy, vol. 39, no. 3–4, pp. 335–346, 2015.
[3] J. Wang, “Charging information collection modeling and analysis of GPRS networks,”
IEEE Trans. Syst. Man Cybern. C Appl. Rev., vol. 37, no. 4, pp. 473–481, 2007
[4] T. Italia, “Telecommunications - SMS, Call, Internet - MI.” Harvard Dataverse, 2015.
[5] G. O’Reilly, A. Jrad, R. Nagarajan, T. Brown, and S. Conrad, “Critical infrastructure
analysis of telecom for natural disasters,” in Networks 2006.pp. 1–6.
[6] Medium. 2021. Decision Trees and Random Forests. [online] Available at: <https://bit.ly/3gHTFqt>

39

JSST2021

Elastodynamic analysis of elastic wave scattering
using coupling method of CQBEM and FEM
Haruhiko Takeda1* and Takahiro Saitoh1
Department of Civil and Environmental Engineering, Gunma University, Japan1*

t160c056@gunma-u.ac.jp

Abstract. The Convolution Quadrature Boundary Element Method (CQBEM) is known as a
new time-domain Boundary Element Method (BEM). The Finite Element Method (FEM) is avery popular simulation tool for various engineering problems. However, the FEM cannot deal
with infinite regions without any modifications. On the other hand, the CQBEM can treat suchregions because the time-domain fundamental solutions used in the CQBEM satisfy the radiation condition automatically. Therefore, in this study, a coupling method of the CQBEM andthe FEM is developed for elastic wave scattering analyses in infinite regions. As numericalexamples, elastic wave scattering in inhomogeneous media is demonstrated using the proposed
method.
Keywords: Finite element method (FEM), Boundary element method(BEM), Time-domain, Convolution quadrature method (CQM),

1. Introduction
The Convolution Quadrature Boundary Element Method (CQBEM) is known as a new timedomain Boundary Element Method (BEM) and effective for wave propagation in homogeneous
regions. The Finite Element Method (FEM) is a very popular simulation tool for various engineering problems and can deal with inhomogeneous regions. However, the FEM generally hastrouble in treating wave propagation in infinite regions. Therefore, it is significant to develop
a new wave propagation simulation tool which can treat both inhomogeneous and infinite regions. Therefore, in this study, a coupling method of the CQBEM and the FEM is developedfor 2-D inplane elastic wave propagation. As numerical examples, elastic wave scattering inhomogeneous media is demonstrated using the proposed method.

2. CQBEM and FEM
As mentione in section 1, the FEM is known as one of the very popular simulation tools. Therefore, the brief description of the CQBEM is discussed in this section. In general, the classicaltime-domain BEM cannot produce stable numerical solutions if we use small time step size.
The CQBEM was developed to overcome this difficulty. The Convolution Quadrature Mehotd(CQM), first proposed by Lubich [1], accurately approximates the convolution integrals by aquadrature formula using the Laplace-transform. The CQBEM is a combined method of the40

JSST2021
time-domain BEM and the CQM. Schantz et al.[2] and Saitoh et al.[3] have solved many wave
problems, such as elastodynamic and viscoelastic wave problems, using the CQBEM.

3. Numerical examples
Some numerical examples are shown in this section. A multiple scattering problem for 2-D
elastodynamics is solved by the proposed method. Numerical analysis model and results canbe seen in Fig.1. 2-D elastic wave propagation in an inhomogeneous and infinite elastic solid
is simulated by the proposed method. The inhomogeneous region inside the dashed square in
Fig.1(a)-(d) is modelized by the image-based method [4] for a concrete block and analyzed by
the voxel FEM. The infinite region outside the inhomogeneous region is treated by the CQBEM.
Two types of materials (aggregates), which have different elastic constants exist in the inhomogeneous region. The continuous boundary condition on the dashed square for the displacement and traction is imposed in order to combine both CQBEM and FEM. As shown in Fig.1,
an incident plane P-wave coming from the infinite horizontal direction propagates the inhomogeneous region, and some of them are scattered and the rest are transmitted. It can be seen thatscattered waves and the incident P-wave propagate to the infinite region without any reflections
because the CQBEM can satisfy the radiation condition.

4. Conclusion
The coupling method of CQBEM and FEM for 2-D elastodynamics was developed in this study.
In the future, the adaptive cross approximation (ACA) method will be applied to save the CPU
time and memory.

Fig.1: Numerical analysis model and 2-D elastic wave fields at various times.

References[1] C. Lubich : Convolution quadrature and discretized operational calculus I, Numer. Math,
Vol.52(1988), pp.129-145.
[2] M. Schanz and H. Antes.: Application of operational quadrature methods in time domainboundary element methods, Meccanica, Vol.32(3) (1997), pp. 179-186.
[3] T. Saitoh, F. Chikazawa and S. Hirose: Convolution quadrature time-domain boundary
element method for 2-D fluid-saturated porous media, Appl. math. Model., Vol.38 (2014),
pp.3724-3740.
[4] K. Terada, T. Miura and N. Kikuchi: Digital image-based modeling applied to the homogenization analysis of composite materials, Comput. Mech., Vol.20(4) (1997), pp.331-346.

41

JSST 2021

Particle-scale FSI computation for internal
ﬂuidization in gravel-particle bed by upward water jet
S. Ushijima1,* , J. Ohno2 , D. Toriu1 , Y. Ueno2
1 Academic Center for Computing and Media Studies (ACCMS), Kyoto University
2 CERE, Graduate School of Engineering, Kyoto University∗ ushijima.satoru.3c@kyoto-u.ac.jp

Abstract. The internal ﬂuidization in a gravel-particle bed by the vertically-upward water jet
entering from the bottom surface was investigated with experiments and the computations
taking account of the particle-scale ﬂuid-solid interactions (FSI). In computations, the ﬂuid
forces acting on about 5,100 gravel particles, each of which is represented with multiple
tetrahedron elements, were calculated from the pressure and viscous terms. As a result, it
was shown that the calculated unsteady characteristics of the internal ﬂuidization and failure
are in good agreement with the experimental results. The consideration was also made interms of the relationships between the ﬂuid dynamic pressure p′ and the eﬀective stress σ′
in the initial shape of the gravel layer.
Keywords: Internal ﬂuidization, ﬂuid-solid interaction (FSI), Parallel computation

1.

Introduction

It has been reported that the serious ground collapses and failures sometimes arise due tothe leakage of the ﬂuids from the underground pipes through which water and oil are transported. Thus, such problems have been investigated from experiments and numerical simulations [1],[2],[3]. In this paper, in order to understand the internal ﬂuidization and resulting
failure of a gravel-particle bed by the vertically-upward water jet entering from the bottom
surface. For that purpose, the unsteady processes are investigated with the experiments andthe particle-scale ﬂuid-solid interaction (FSI) numerical prediction using parallel computations. The predicted results were compared with the experimentally obtained behaviors. The
mechanisms of the internal ﬂuidization was discussed using the numerically obtained ﬂuid
forces acting on the particles and pore ﬂuid pressures.

2.

Outline of Experiments

Figure 1 shows the outline of the main part of the experimental setup and computational
area. The lengths indicated in Fig. 1 are as follows: l1 = 350, l2 = 40 and l3 = 200 [mm],42while hg ≈ 100 [mm] and hw ≈ 300 [mm]. The inner diameter of the bottom pipe D is 31[mm], from which the water jet of the average velocity 0.21 [m/s] enters in the verticallyupward. As shown in Fig. 1, gravel particles of the average diameter 7 [mm] were ﬁlled in
the gravel box. The approximate depth of the gravel bed is 100 [mm] and the porosity is
about 0.44.
computational area

l2

l1
x3

gravelbox
gravelparticles hg

x1

hw
l3

D
flow

x3
x2
flow

Fig.1: Experimental tank and computational area (left = front-, right = side-view)

3.

Particle-scale FSI computations

The ﬂuidized processes of the gravel bed were numerically simulated using the computation
method proposed in [4]. Figure 2 (a) illustrates the Eulerian cells used for ﬂuid computations
as well as the Lagrangian unstructured cells representing a gravel particle. As schematicallyshown in Fig. 2, the Eulerian cell size is 0.6 [mm], which is suﬃciently ﬁner than the gravel
particle size (7 [mm]), so that the surrounding ﬂows around a particle can be resolved.
The total number of the gravel particle is about 5,100, each of which is represented with
about 100 tetrahedron elements. The massive computations were executed with 272 parallel
computations with the domain decomposition as shown in Fig. 2 (b).

(a) Eulerian and Lagrangian cells(b) Domain decomposition for parallel computation

Fig.2: Computation cells and parallel computation
Figure 3 shows the comparisons between experiments and computations for the unsteadyprocesses from initial condition to the ﬂuidization of gravel particles and ﬁnal penetrationof the jet ﬂow throughout the bed. The color contours in the calculated results show themagnitude ω of the vorticity vector of the ﬂuid. The range of ω shown in Fig. 3 is between
0 [1/s] (blue) and 600 [1/s] (red).

43

(a) experiment: t = 0.00 [s]，computation: t = 0.00 [s](b) experiment: t = 0.51 [s]，computation: t = 0.495 [s](c) experiment: t = 1.50 [s]，computation: t = 1.50 [s](d) experiment: t = 3.02 [s]，computation: t = 3.00 [s]Fig.3: Comparisons of gravel particles (left = experiments, right = computations)Figure 4 shows the distributions of positive p′ − σ′ , where p′ is pore dynamic pressure
among particles and σ′ is the eﬀective stress in the initial gravel bed shape calculated withthe speciﬁc gravity (about 2.59) and the porosity (about 0.43). The range of p′ − σ′ shownin Fig. 4 is between 0 [Pa] (blue) and 600 [Pa] (red). Figure 4 also shows the ﬂuid forces
F F = |F| calculated from the pressure and viscous terms in the ﬂuid computations. It can beseen that the ﬂuidized areas, where p′ − σ′ and F F are relatively high, are closely related tothe unsteady behaviors of gravel particles in the bed.

44

(a) p′ − σ′ (t = 0.045 [s], left) and F f (t = 0.15 [s], right)(b) p′ − σ′ (t = 0.495 [s], left) and F f (t = 0.495 [s], right)(c) p′ − σ′ (t = 1.50 [s], left) and F f (t = 1.50 [s], right)
Fig.4: Distributions of p′ − σ′ (left) and F F (right)

References[1] Alsaydalani, M. and Clayton, C.: Internal ﬂuidization in granular soils, Journal of Geotechnicaland Geoenvironmental Engineering, Vol. 140, No. 3, p. 04013024, 2014.
[2] He, Y., Zhu, D. Z., Zhang, T., Shao, Y. and Yu, T.: Experimental observations on the initiationof sand-bed erosion by an upward water jet, Journal of Hydraulic Engineering, Vol. 143, No. 7,
p. 06017007, 2017.
[3] Ngoma, J., Philippe, P., Bonelli, S., Radjaı̈, F. and Delenne, J.-Y.: Two-dimensional numericalsimulation of chimney ﬂuidization in a granular medium using a combination of discrete element
and lattice Boltzmann methods, Physical Review E, Vol. 97, No. 5, p. 052902, 2018.
[4] Ushijima, S., Toriu, D., Yanagi, H. and Tanaka, H.: Numerical prediction for transportation ofgravel particles and saltation-collapse equilibrium due to vertical jet, Journal of Applied Mechanics, JSCE Ser.A2, Vol. 75, No. 2, pp. I 289–I 300, 2019.

45

JSST2021

Embodiment Design of Biped Robot for Walking on
Rough GroundKakeru Yokoyama1*, Van-Tinh Nguyen2, Hiroshi Hasegawa1
Systems Engineering and Science, Graduate School of Engineering and Science,

1

Shibaura Institute of Technology, Japan
School of Mechanical Engineering, Hanoi University of Science and Technology, Vietnam*mf20083@shibaura-it.ac.jp

2

Abstract. Biped robots consume a lot of energy due to the attitude control during walking,which causes the running time to be shortened. The robot that simulates arm swing walking,which is a passive human movement, are optimized by simulating with mechanical analysis
software without considering joints and actuators. Therefore, the design on details of mechanism of the biped robot is expected for this work. As a result of simulation, the model of thisstudy generally walked stably, however it is necessary to improve the width of the body in thefront-back direction for further walking stability.
Keywords: Embodiment design, Multibody dynamics analysis, Biped robot

1. Introduction
Many biped robots have been researched and developed to support human daily life in a moresophisticated way, and some of them can walk with a stable behavior even on rough ground.
On the other hand, their energy consumption is high, which is one of the reasons for the decrease in running time, since actuators such as servomotors are frequently used for posturalcontrol. To solve this problem, V.T. Nguyen et al. focused on the passive motion of the humanupper body to optimize the robot's gait, and developed a basic model in which the shoulder andthe opposite hip joint are connected by a spring damper [1]. This basic model was designed
without considering the detailed mechanisms of the actuators and joints and the interference ofthe mechanisms, and it is difficult to reflect the model directly in actual machines. In this study,
we will perform an embodiment design that can be fabricated from the basic model [1].

2. Walking simulation for embodiment design
The basic model [1] and the designed embodiment model are shown in Figure 1. A multi-dynamics simulation is performed to confirm the walking behavior of the embodiment model.
For this simulation, we use the mechanical analysis software Adams from MSC Software. For
the walking behavior, a gait function derived from the human walking pattern is used. The
characteristics are expressed and formulated by using triangular functions from the data of the46

JSST2021
human walking cycle [2]. The gait function is shown in Equation (1).
φ (𝑡) = 𝑎 + 𝑏 cos(𝜔𝑡) + 𝑐 sin(𝜔𝑡) + 𝑑 𝑐𝑜𝑠(2𝜔𝑡)(1)where φ is the angle of the joint, t is the time, ω is the angular velocity, and ai, bi, ci, and di
are the coefficients.
By changing these coefficients, we can derive the gait function to be assigned to the embodiment model. In this study, we use the same parameters as in the basic model [1]. The angle
trajectories of each joint in this study are shown in Figure 2.

Figure 1 Right: basic model. Left: embodiment model

Figure 2 Angular trajectories

47

JSST2021

3. Results and Discussion
Figure 3 shows the graphs of the lateral and vertical displacements of the head in the gait of
the basic model [1] and the embodiment model. The horizontal direction shows that the difference in displacement between the basic model [1] and the embodiment model is larger. This isdue to the fact that the embodiment model is wider in the front-back direction than the basicmodel [1], which tends to swing back and forth during walking. It is considered that this resulted in an unstable posture on the rough ground, which caused the robot to move in a slightly
slanted direction. On the other hand, in the vertical direction, the period of the displacement
was constant after 4 seconds, which is the period of the flat ground, and it can be concluded
that the robot of the embodiment model walked stably. Using the information obtained fromthis study, we believe it is necessary to make improvements in order to reduce the width of therobot in the front-back direction additionally.

(a) lateral direction(b) vertical direction
Figure 3 Comparison of head displacement during walking

References[1] V.T. Nguyen, D. Kiuchi, H. Hasegawa: Development of Foot Structure for Humanoid
Robot Using Topology Optimization, Tech Publications, Advanced Engineering Forum,
volume:29 (2018).
[2] N. Ito, H. Hasegawa: The Robust Design to Generate the Gait Pattern of a Small Biped
Robot, Japan Society for Design Engineering, Design Engineering, volume: 45, (6),(2010) (in Japanese).

48

JSST 2021

Effect of Adsorption and Reactive Behavior of Water
Vapor in Reactive Magnetron Sputtering Simulations
using Berg’s Model
Allen Vincent B. Catapang1,* and Motoi Wada1
1 Graduate School of Science and Engineering, Doshisha University∗ cyjf3301@mail4.doshisha.ac.jp

Abstract.
The time-resolved Berg’s model was developed as a first-order model for describing thesurface processes occuring at different regions in a reactive magnetron sputtering discharge.
In this study, the effect of using a highly adsorptive gas, water vapor, as reactive gas inBerg’s model is investigated. The contributions due to parameters affecting the adsorption
onto the surface, the flow rate, Qr , sputtering yield ratio of metal to compound, YYmc , sticking
coefficient, α, and surface site density, n s , were investigated and the results were comparedto obtained experimental data.
Keywords: Reactive Magnetron Sputtering, Water Vapor, Berg’s Model

1.

Introduction

As compared to sputtering with inert gases, the presence of a reactive gas in plasma often
results to complex behavior, such as the hysteresis of process parameters and discharge
instabilities. To explain the nonlinear plasma behavior, multiple models have been explored,such as Berg’s model, which has served as a basic, fundamental model, and the Reactive
Sputter Deposition (RSD) model [1, 2].
Water vapor has been used in plasma-based processes, such as surface modification andthin film deposition. In deposition processes, it is used as a substitute to conventional oxygengas as well as a measure to introduce hydrogen into the plasma. Particularly in semiconducting oxides, the added hydrogen acts as a shallow donor dopant, improving film conductivity
and reducing crystal lattice stress without extensively changing the band structure of the film[3]. However, the behavior of water vapor in vacuum processes can be complex, since there
can be multiple interactions involved [4]. Modelling of the process would clarify and explain the observed phenomena involved, allowing for a more widespread use of water vapor
in film deposition.
Shown in Fig. 1 is the time-resolved pressure as the plasma is ignited, and fitting ofBerg’s model to experimental data. An initial drop in the pressure is observed, followed bya gradual, successive recovery towards a saturation level. In this paper, Berg’s model wasused to investigate the water vapor plasma and explain the observed change in pressure, with
an emphasis on the parameters that are highly affected by reactive and enhanced adsorption
behavior.

49

Figure 1: Measured pressure of the reactive magnetron sputtering discharge at varying initialpressures of pure water vapor, at 100 mA discharge current and a 70mm-φ metallic Zn(99.2%) target, and the fitting of Berg’s model to the pressure at 1.0Pa.

2.

Simulation ModelIn this study, the model used is Berg’s time-resolved model, one of the fundamental, firstorder models used in reactive magnetron sputtering [2]. This model simplifies some of thecomplicated but important phenomena, however, it has been able to identify key processing problems, such as the hysteresis of discharge parameters [1]. Thus, this model can bereasonably used to identify the key mechanisms involved when water vapor is utilized inreactive magnetron sputtering.

Figure 2: The processes involved in the time-resolved Berg’s model
As shown in Fig. 2, the pressure of a reactive magnetron sputtering system can be
described by the flow of gas in the reactive sputtering chamber, where Q is the gas flow rate,with the subscript denoting the region (r: reactive gas flow, s: substrate, t: target), v is thechamber volume (in m3 ), kb is Boltzmann’s constant, and T is the temperature. For Q pump ,50

3

the flow rate depends on the pumping speed, S (in ms ).
p(t)Fr (T ) = √2πmr kb T(1)

The values of Q s and Qt are determined by the flux towards the surface, Fr , obtained
from the kinetic theory of gases in Eq. 1, the sticking probability of the reactive gas moleculeon the surface, α, the area, A (in m2 ), and the compound fraction at the surface, θ. The valueof θ is important, as the presence of the compound determines the gettering of the reactivegas molecule as it interacts with the metallic particle. At the target surface, this describes the
poisoning phenomena, while at the substrate, it determines the uniform film composition. In
Fig. 2, the value of θ is dependent on the surface site density, n s , with a base value related tothe atomic surface density of the material, the flux towards the surface, and the removal or
deposition of material, determined by Jion , the ion current density, and the sputtering yield,
Y. The parameter z is the stoichiometric factor of the oxide compound. By solving theequations simultaneously, the change in pressure for reactive magnetron sputtering can be
determined [1, 2].

3.

Results

To understand the effect of the parameters, the difference of water vapor to conventional
gases must be established. First, the dissociation of the water molecule can result to distinct
species, O, H and OH. These molecules behave differently in sputtering and film formation.
Second, water easily adsorbs onto a surface. Multilayer adsorption usually occurs, with adesorption energy ranging from 0.8 to 8 eV, depending on the type of bond and the thicknessof the adlayer.
Shown in Fig. 3 are the effects of specific parameters determined to be affected by thedifference in behavior. First is the flow rate of the reactive gas. As the flow rate is increased,
the initial drop in pressure becomes more prominent, since there is an increased flux ofmolecules that can contribute to the surface processes that consume the reactive gas.
The second parameter is the ratio of the sputtering yield of the metal and the compound. At decreased values of the compound sputtering yield, other surface process can be
accounted for, such as redeposition or readsorption. At low ratios, the saturation pressure
increases due to a decreased contribution of Q s and Qt as the compound is formed at thesurfaces, indicating decreased gettering when more compound is present.
The third parameter is the sticking coefficient, α. The magnitude of the coefficient determines how tightly a molecule will adhere on the surface, and given the concentration ofplasma near the target, the value of αt is expected to be less than that of α s due to energetic
particle bombardment. From the fitting of the experimental data at 1.0Pa, αt was at 9x10−5 ,while α s is at 1x10−4 . As shown in Fig. 3, higher values of αt results to a more prominent
drop in the pressure, since it influences the consumption of gas at the target, while higher
values of α s increases the saturation pressure level, due to increased compound formation at
the substrate area.
Lastly, the surface site density, n s is changed to account for the possible multilayer
adsorption of the molecule at the surface, in the form of k(n s ). At the initial pressure considered, water vapor favors the formation of multiple layers through physisorption, leading toan increased site density for compound formation. As seen in Fig. 3, the time at which the51

Figure 3: Effect of the parameters, flow rate (top left, 0.11-0.27sccm), YYmc (top right, 0.0010.009), α (bottom right, 2x10−5 -1x10−4 ), and n s (bottom left, 400-1000 layers) in the timeresolved Berg’s model. The default values for the parameters are: flow rate=0.19sccm,
Yc−5−4
Ym =0.003, αt = 9x10 , α s =1x10 , and n s = 850 layers.
pressure minima occurs shifts and the rate of increase to the saturation pressure is delayed
at a higher number of layers, since there are more reaction sites for compound formation.
By combining the effect of the parameters, the observed pressure of water vapor plasma
can be explained, and an appropriate fitting of the experimental data to the model can be
performed. By using this model, it allows us to predict the behavior of the plasma, and thefilm formation process when water vapor is used.

References[1] K. Strijckmans, R. Schelfhout, D. Depla: Tutorial: Hysteresis during reactive magnetron sputtering processes, Journal of Applied Physics, 124 (2018), 241101-1–
241101-30.
[2] S. Berg, T. Nyberg: Fundamental understanding and modeling of reactive sputtering
processes, Thin Solid Films, 476 (2005), 215–230.
[3] C.G. Van De Walle: Hydrogen as a cause of doping in zinc oxide, Physical Review
Letters, 85 (2000), 1012–1015.
[4] A. Berman: Water vapor in vacuum systems, Vacuum, 47 (1996), 327–332.

52

JSST 2021

Target ionization and surface damage models ofgraphite via nanosecond pulse irradiation
James Edward Hernandez1,* , Kengo Moribayashi2 , James Koga2 , Motoi Wada1
1 Graduate School of Science and Engineering, Doshisha University
2 Kansai Photon Science Institute, National Institutes for Quantum and Radiological

Science and Technology∗ cyjd3302@mail4.doshisha.ac.jp

Abstract. A numerical model of nanosecond pulsed laser ablation on a graphite surface isbeing developed. The rates of photoionization of carbon atoms released by the ablation areevaluated. Total ion energy of irradiated atoms are calculated via photoionization rates using
Monte Carlo simulation.
Keywords: Nanosecond laser, Laser ablation, Numerical Modelling

1.

Introduction

Nanosecond laser ablation is widely employed in pulsed laser deposition and mass spectrometry [1, 2]. Ablation of the target via a nanosecond laser involves several mechanisms
such as ion formation, target heating, and particle expulsion, which are heavily dependent
on the laser parameters. In the initial phase of laser ablation, primary interaction involves
photoionization of target bound states via the laser field, followed by plasma formation andcollisions between electrons and atoms [3]. In nanosecond laser ablation, the electron emission is described to be thermionic, rather than dominated by photoionization [4]. Following
the electron production, ions are formed from the surface due to electron-atom collisions,and a plasma is formed which is assumed to follow a Maxwellian distribution [5]. The
plasma formation implies that the temperature exceeds the vaporization threshold of the target. Simulations are performed to quantify the target ablation rate evolution using thermal
model, whereby the ablation rate stabilizes after tens of ns following the laser pulse [6]. The
aforementioned models refer to the case where the target is at a previously unablated state(i.e. a fresh target). In this work, electron production in the geometry of a formed cavity isinvestigated, focusing on the plasma formation in the context of ion and electron generation.
Electron production is simulated using Monte Carlo methods [7]. The simulation employsa vacuum environment.

53

2.

Numerical model

The simulation is set up in 3D coordinates with x is specified along the target surface normal.
A 1064 nm wavelength laser with full width half maximum (FWHM) of 5 ns and 5 GW/cm2
laser intensity is incident on a graphite target at x = 0, y = 0, z = 0. The pulse peak islocated at t = 30 ns. The timestep used is 1 ps. Atoms are placed in 3D coordinates eachwith spacing dr = 0.2 nm to simulate graphite density, with n2 atoms placed in the x-planeand n1 × n1 atoms in the y-z plane. Here n2 = 10, n1 = 20 is employed. The atoms arelocated on the target surface are arranged in the y-z plane labelled from 1-100, and 101-200
for x = dr, and so on.
The incident laser is represented by a Gaussian source term and is described as
S (x, y, z, t) = Io e(−x/λD ) −e(t−t peak )22σ2te− (y−yc2)
2σy

2

2

e

− (z−zc2)2σz(1)

where Io is the peak laser intensity, λD = 1.5 dr is the decay rate of the laser intensity along
the target surface normal. First, we show the ionization probability due to multiphoton
ionization (MPI). The probability for an electron to change the level from a ground state toan excited state is given by [8]5/2!3/2 
!
B + ΦpBγ S γ,P MPI = ω
 p
Eω
Eω
1 + γ2p 2(B + Φ p ) 
1 + γ2 
−1
exp −sinh γ − γ ,Eω
1 + 2γ2 


∞
X


γ
−1
×
S (γ, x) =exp −2[hx + 1i − x + n] sinh γ − p
2
1
+
γ
n=0
v
t

  2γΦ   p(hx + 1i − x + 1)
1 + γ2(2)pwhere Eω = ~ω is the photon energy, γ = B/(2Φ p ) is the Keldysh parameter, with B, Φ p =e2 E 2 /4mω2 being the ionization potential and potential energy. Angle brackets represent theconversion to integer value. The function Φ(x) is the Dawson integral. [9]The laser strikes the target with a FWHM defined by σy and σz in Eq. 1 corresponding toa diameter of 2 nm. The penetration of the laser depends on λD assumed to be 1.5 dr. Whenthe laser irradiates the target, the atoms ionize according to P MPI . A random number RN isgenerated for each timesteppsuch that when RN < P MPI the atom ionizes and a photoelectronis ejected with a velocity 2(nhν − φ)/me with a random direction.√Electrons are assumed
to follow Maxwell-Boltzman characteristics, so the velocity is v = kb T/me . Electrons areaccelerated due to the laser field via the ponderomotive force whose magnitude is a pond =e2 m2e E 2 /ω2 , where e,me , and ω are the electronic charge, mass, and laser frequency, andlaser electric field respectively. Electron-electron collisions reduce the particle velocity byme~vνee , where νee is the electron-electron collision frequency.

54

Fig. 1 shows the temporal laser pulse shape with the P MPI for varying laser intensities.
The power law dependence of the P MPI with the laser intensity is shown.

Figure 1: (left) Laser intensity and P MPI for varying laser intensity, (right) Dependence ofP MPI on laser intensity

Figure 2: (left) Time evolution of ionized atoms for varying laser intensities, (right) maximum number of ionized atoms versus peak laser intensity

3.

Results

Figure 2 shows the time evolution of the maximum number of photoionized atoms ni,max forvarying laser intensities. Increase in ion number is observed with laser intensity, indicating
the higher number of available photons for ionization. As the laser intensity increases, thesaturation of the ion number occurs earlier relative to the laser pulse. This number is plottedagainst the peak laser intensity, as shown in the right figure of Fig. 2. Rapid increase inionization rate is observed for lower laser intensities. Fig. 3 shows the atom label vs timefor varying laser intensity. Large number of atoms are ionized on the surface, and at thecenter (vicinity of n1 = 50) even before the time for the laser to reach its peak intensity. Thisshows the preferential direction of ionization at the target center. Figure 4 (right) shows thetotal target energy as a function of time. Ion energies rapidly increase at an earlier phase forincreasing laser intensity. Future work will incorporate motion of the ionized atoms which
is neglected in the current model.

55

Figure 3: (left) Atom label vs time for varying laser intensity, (right) target total energy vstime for varying laser intensity

References[1] H. Jadhav, A.K. Singh, S. Sinha: Pulsed laser deposition of graphite in air and in
vacuum for field emission studies, Mat. Chem. Phys., 162 (2015), 279-285.
[2] A. O’Keefe, M. M. Ross, A. P. Baronavski: Production of large carbon cluster ions bylaser vaporization, Chem. Phys. Lett., 130:1, 2 (1986), 17-19.
[3] S. Amoruso, R. Bruzzese, N. Spinelli, R. Velotta: Characterization of laser-ablation
plasmas, J. Phys. B: At. Mol. Opt. Phys., 32 (1999), 131-172.
[4] J. F. Ready: Effects of High Power Laser Radiation, Academic Press, (1971), 134-135.
[5] M. Aden, E. W. Kreutz, A. Voss: Laser-induced plasma formation during pulsed laser
deposition, J. Phys. D: Appl. Phys., 26 (1993), 1545-1553.
[6] S. Sinha: Nanosecond laser ablation of graphite: A thermal model based simulation,
J. Laser Appl., 30 (2018), 012008.
[7] K. Moribayashi: Spherically symmetric models for x-ray damage and the movementof electrons produced in non-spherically symmetric targets such as bio-molecules, J.
Phys. B: At. Mol. Opt. Phys., 43 (2010), 165602.
[8] Y. Kishimoto, S. Kato: Modeling of Ionization in Particle Simulation, J. Plasma
Fusion Res., 90, 7 (2014), 412-418.
[9] F. G. Lether: Constrained near-minimax rational approximations to Dawson’s integral,Appl. Math. Comput., 28, 2-3 (1997), 267-274.

56

JSST 2021

Analysis of the space charge effect in an ion mobility
spectrometer
Keith Nealson M. Penado1,* , James Edward A. Hernandez1 , Allen Vincent B.
Catapang1 , Motoi Wada1
1 Graduate School of Science and Engineering, Doshisha University∗ cyjf3303@mail4.doshisha.ac.jp

Abstract. Ion mobility spectrometry has been a useful tool in the detection of dangerous
substances such as explosives and narcotics. The application of this method has expanded
to the analysis of atmospheric pressure plasma. Under ambient conditions, the density ofgas particles is much greater than the density of ions present in the spectrometer such that
Coulomb’s forces do not dominate as much. However, the space charge contributions of
the ions to the mass analysis electric field should affect the resolution of the spectrometersignal. This study aims to analyze this effect by using a collisional model based upon thehard-sphere model and the Monte Carlo method. The results of the simulation will then becompared with experimental results to determine the accuracy of the model.
Keywords: Plasma physics, Monte Carlo method, Computer simulation

1.

Introduction

Ion Mobility Spectrometry (IMS) is an electrophoretic analysis method used to identify ionsin a swarm based on their mobility, K, in a buffer gas. The mobility is related to the drift
velocity, vd , of the ion as it moves through the spectrometer guided by the electric field, E[1];
vd = KE(1)

The intrinsic characteristics of the ions such as mass, charge, and structure affect its
mobility. The combination of such characteristics interacting with the electric field and gas
results to the separation of the swarm into discrete groups which would be analyzed at theend of the spectrometer [2]. If the density of ions that enter the spectrometer is much lessthan the density of the gas particles, then the Coulomb interactions between the ions will notdominate. However, the high density ions should affect the local space charge of the system
which could influence the drift behavior of the swarm as they move across the IMS. Thisphenomenon will be clarified by simulating the behavior of Ar+ ions in the IMS using thehard-sphere collisional model and the Monte-Carlo method.

57

Figure 1: Schematic drawing of the IMS device (left) and the corresponding potential field
distribution (right) along the z = 0 plane. The IMS device shows (a) the opening orifice, (b)the electrode plate, (c) the insular spacers, (d) the end plate detector, and (e) the cyclindricalshielding. The potential distribution was calculated using the AMaze software.

2.

Computational Method

The IMS device used in this study exhibits a local electric field of E = 300 V/cm at thecenter under ambient temperature and pressure (T = 273.15 K, P = 1 atm) with air being
the buffer gas. The collisions between ions and gas particles are treated as hard spherecollisions [3, 4], with the probability, p, being a reduced exponential function of the ratio ofthe distance travelled by the ion, d f , and its respective mean free path, λ (Eq.2). The mean
free path, λ, is given by a standard equation (Eq.3) with pressure, P, temperature, T , and,the radius of the gas particle, d, as its factors. The simulation assumes that a collision occurs
if a random number between 0 and 1 is greater than the probability, p.
p = 1 − exp (−d f /λ)

(2)

kb T
(3)
√
πd2 P 2The current simulation utilizes Coulomb’s law to approximate the space charge contribution of an ion by defining a voltage threshold of 2 µV from an individual point charge. To
further clarify the space charge effect, the current model multiplies the voltage contributionof each particle by a factor of 10. This voltage contribution is added to the background
potential matrix before allowing the particles to move for each time step iteration.
The ion trajectories and time of flight signal will be obtained for the case with the inclusion and absence of the space charge effect. The mobilities of the ions will also be comparedbased upon the results of the simulation.
λ=

58

3.

Results and Discussion

The ion trajectories (Fig.2) and time of flight signal (Fig.3) for the movement of particles in
the absence, and presence of the space charge contributions for Ar+ ions are shown in the
figures below. The time steps used in these cases were chosen such that any further decreasein the time step will not result to any change in the peak drift time.

Figure 2: Ion trajectories of 10 Ar+ ions with (left) and without (right) space charge contributions. In the absence of a space charge, ions were made to interact only with the electricfield such that their trajectories are independent from each other.
The number of particles detected at the end of the spectrometer is greatly reduced withthe addition of the space charge effect. Peak reduced mobility values for both cases however,is around 1935 cm2 /Vs. The resolution of the spectra varies; with the FWHM values being
4.98 µs, and 1.05 µs for the spectrum with and without charge contributions respectively.
This suggests that the space charge effect of ions results to the broadening of the spectral
peak. Increasing gas temperature from 273 K to 300 K in the case with the space charge
contributions results to a sharper spectrum. The FWHM value at a gas temperature of 300
K is around 4.78 µs.

Figure 3: Time of flight signal of 10,000 Ar+ ions with and without space charge contributions at a spectral resolution of 0.15 µs. Peak time of flight for both cases is around 51.67µs.

59

Figure 4: Time of flight signal of 10,000 Ar+ ions with the space charge contribution at gas
temperatures of 273 K, and 300 K. Spectral resolution is 0.15 µs

Further investigation with other models to approximate the space charge will be needed
to investigate the actual behavior of particles in the experimental system since the current
simulation could only compute up to 104 particles efficiently. With an increase in the numberof particles, the multiplicative modification to the space charge contribution can be removed
which will result to a more accurate representation of the system. The results of the simulation will also be compared to actual experimental results to determine the validity of themodel used in this study.

References[1] G. Eicemann, Z. Karpas, H. Hill Jr.: Ion Mobility Spectrometery (3rd Ed.), Parkway
NW: CRC Press (2014).
[2] S. Harvey, R. MacPhee, C. Barran: Ion mobility mass spectrometry for peptide analysis, Methods, 54 (2011), 456–461.
[3] J. Xu, W. Whitten: Monte Carlo simulation of ion transport in ion mobility spectrometry, Int. J. Ion Mobil. Spec., 11 (2008), 13–17.
[4] A. Appelhans, D. Dahl: Measurement of external ion injection and trapping efficiency
in the ion trap mass spectrometer and comparison with a predictive model, Int. J. Mass
Spec., 216 (2002), 269–284.

60

JSST 2021

脱塩基部位からなるクラスターDNA 損傷の分子動力学シミュレーション
Molecular Dynamics Simulation of Clustered DNA
Damage Composed of Apurinic/Apyrimidinic Sites
寺川 和志 1*, 川波 竜太 1, LI HAOLUN1, 藤原 進 1, 水口 朋子 1, 米谷 佳晃 2,
中村 浩章 3,4Kazushi Terakawa1*, Kawanami Ryuta1, Haolun Li1, Susumu Fujiwara1, Tomoko Mizuguchi1, Yoshiteru Yonetani2, Hiroaki Nakamura3,4

1

京都工芸繊維大学,
2

QST,

3

名古屋大学,
4

核融合研

Kyoto Inst. of Tech.,

1

2
3

QST,

Nagoya Univ.,
4

NIFS

*

m1672015@edu.kit.ac.jp

Abstract. DNA はアルキル化や酸化、放射線照射などに弱く、これらを受けると損傷してしまう。DNA 損傷には、8-オキソグアニン(8oxoG)、脱塩基部位(AP 部位)、一本鎖切断（SSB）、二本鎖切断（DSB）などがある。複数の DNA 損傷が局所的に発生（クラスター化）すると、細胞中の酵素による修復が困難となるため細胞死やがん化の原因になる。AP 部位は DNA 中の塩基が脱離している状態であり、DNA に大きな影響を与える SSB や DSB につながる。本研究では、AP 部位からなるクラスター
DNA 損傷の構造変化を明らかにするため、AP 部位を有する損傷 DNA の分子動力学シミュレーションを行い、DNA の構造変化を解析した。
Keywords: molecular dynamics simulation, apurinic/apyrimidinic (AP) site, clustered DNA
damage61

JSST 2021

1. Introduction
DNA はデオキシリボースとリン酸からなる主鎖と、リン酸に結合している塩基から構成される多くの生物の遺伝情報の継承と発現を担う物質であり、２本のポリヌク
レオチド鎖は相補的な塩基対(A-T 塩基対、G-C 塩基対)間の水素結合により結合し、右巻き二重らせん構造をとっている。
DNA はアルキル化や酸化、放射線照射などに弱く、これらを受けると損傷してしまう。特に 2011 年の福島原発事故で、放射線汚染が問題となった。放射線は細胞に当たると、細胞の中にある遺伝子の本体である DNA に傷をつけることがある。DNA
損傷には DNA のグアニン塩基が酸化されてしまう 8-オキソグアニン(8oxoG)、塩基が脱離してしまう脱塩基部位(AP 部位)、DNA 二重らせん構造の片方が切れてしまう
一本鎖切断(SSB)または両方が切れてしまう二本鎖切断(DSB)などがある。これらの損
傷が局所的に集まることをクラスター化と呼ぶ。これらの傷は通常、体内の酵素により修復されるが、傷が多かったり、傷がクラスター化すると酵素による修復が不十分になり、細胞死やがんに繋がることがある。DNA の相補鎖内でのクラスターDNA 損
傷の分子動力学(MD)シミュレーションは行われており[1]、本研究では同一鎖内でのクラスターDNA 損傷について MD シミュレーションを行い、DNA の構造の変化を解析した。

2. Models & Simulation
以下のように未損傷 DNA と損傷 DNA に対して平衡化を行った。シミュレーショ
ンソフトは分子設計、MD 計算、計算結果の解析などを行うことができる、Amber20を用いた。
2.1 未損傷 DNA
Amber20 の Nucleic Acid Builder モジュールを用いて未損傷の DNA を作成した。
次に、Amber20 の xLEAP モジュールを用いて、電荷を中和するために対イオン (Na+)
を入れ、DNA の周りを溶媒(水)で満たした。DNA の力場は Amber20 の bsc1 力場を、
水の力場は TIP3P 力場を用いた。続いて、溶媒と対イオンのエネルギーを最小化した
後、系全体のエネルギーを最小化した。得られた配置を初期配置として、一定体積の周期境界条件の下、系の温度を 0 K から 310 K まで昇温した。(20ps, 10,000 ステッ
プ, 1 ステップ=2 fs) 。そして、温度を 310 K に維持し、一定圧力(1 atm)の周期境界条件の下、9ns (4,500,000 ステップ)の平衡化を行った。
2.2 損傷 DNA
このようにして平衡化した未損傷 DNA(Fig.1 (a)) に、脱塩基部位を１か所または２か所挿入した(Fig.1 (b))。Amber20 の xLEAP モジュールを用いて、電荷を中和するために対イオン (Na+)を入れ、DNA の周りを溶媒(水)で満たした。DNA の力場は AP 部
位に関する補正を含む bsc1 力場を、水の力場は TIP3P 力場を採用した。対イオンと溶媒を加えてそれらのエネルギー最小化を行い、その後系全体のエネルギーを最小化した。その配置を初期配置として、一定体積の周期境界条件の下、系の温度を 0K から 310K まで昇温した(100ps, 50,000 ステップ)。さらに、温度を 310K に維持し、一定圧力(1atm)の周期境界条件の下、25ns の平衡化を行った (12,500,000 ステップ)。上
記のプロセスを Fig.2 に示す 7 種類の損傷 DNA に対して行った。損傷 DNA の構
造変化を調べるため、損傷 DNA の平均二乗偏差(RMSD)と溶媒接触表面積(SASA)の解析を行った。また、比較のため未損傷 DNA に対しても同様の解析を行った。

62

JSST 2021
(a)

(b)

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

Fig.2 Models of clustered DNA damage composedof 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.

3. Results
3.1 Root Mean Square Deviation
AP 部位による構造の変化を解析するため、平均二乗偏差(RMSD)を計算した。
RMSD(Root Mean Square Deviation)は以下の式で表され、ある構造｛𝒓𝑖 ｝と参照構造｛𝒓0𝑖 ｝とのずれを記述する。構造がどれだけ変化したかの指標になる。RMSD の値が大きいとより構造が変化しているといえる。𝑖は計算した原子の番号を表し、𝑁は計算した原子の個数を表す。
𝑁

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

RMSD [Å]DNA 主鎖（P, O, C を対象）についての RMSD を 2ps ごとに計算した。未損傷 DNA
及び損傷 DNA の時間に対する RMSD を Fig.3 に示す。

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

63

JSST 2021
Fig.3 によると、未損傷 DNA の RMSD(undamaged)より、損傷 DNA の RMSD(seq0,
seq±1, 3, 5)が大きい傾向にあると判断できる。全時間の RMSD の平均値の大きさを比較すると、seq0 < undamaged < seq-3 < seq+3 < seq +5 < seq+1 < seq-1 < seq-5 となった。
seq-5 では主鎖の曲がりが他の DNA より大きくなっているのを確認でき、RMSD の平均値が大きくなっている原因だと推測される。損傷 DNA は未損傷 DNA より構造が変化しているといえる。
3.2 Solvent-Accesible Surface Area
OH ラジカルが DNA に接近することで、水素引き抜きなどにより DNA に損傷が生じる。OH ラジカルの電気双極子モーメントは H2O の電気双極子モーメントと値が近く静電特性が似ている。また、水の接近数と溶媒接触表面積(SASA)は相関があるた
め、H2O の SASA を計算することによって OH ラジカルによって損傷を受けやすい部位の予測をすることができる[2]。SASA(Solvent-Accessible Surface Area)とは、溶質分子が溶媒に接触する曲面の面積のことであり、SASA の値が大きいほど OH ラジカ
ルとの接近回数が多くなり、損傷を受けやすいと予測される。
今回は、DNA 主鎖の H1’, H2’1, H2’2, H3’, H4’, H5’1, H5’2 の 7 つの H 原子の SASAの時間平均を計算した。H2’=H2’1+H2’2, H5’=H5’1+H5’2 とした。H 原子からの距離を 1.4Å に指定して計算を行った。各 DNA の SASA を Fig.4 に示す。
Fig.4 によると、損傷 DNA の SASA は未損傷 DNA と比べて、塩基脱離した糖の H2’
は値がかなり大きくなっており、H4’と H5’は値が小さくなっていた。したがって、塩基が脱離している部分は溶媒と接触しやすいため、損傷を受けやすいと予測できる。

Fig.4 The values of SASA for an undamaged DNA and various damaged DNAs. The verticalaxis represents SASA in units ofÅ2. The upper panels represent strand1 and the lower panels
represent strand2 for each DNA.

4. Future Work
引き続き MD 計算を行い、DNA の構造変化に対して、その変化を特徴づける物理量を定義し、定量的に評価を行う。また、構造の変化、SASA の変化と糖の C1’間距離に相関があるかどうか評価を行う。

References[1]

E. Bignon, H. Gattuso, C. Morell, F. Dehez, A. G. Georgakilas, A. Monari, and E. Dumont(2016). Nucleic Acids Research. 44, 8588-8599[2] Y. Yonetani, H. Nakagawa (2020). Chemical Physics Letters. 749, 137441

64

JSST2021

環境を考慮した中国での製造に関する検討
Consideration of environment-friendly
manufacturing in China
周沛涵, 杉本 等，大塚晃

ZHOU PEIHAN, Hitoshi Sugimoto, Akira Otsuka
事業創造大学院大学 〒950-0916 新潟県新潟市中央区米山 3 丁目 1-46
Graduate Institute for Entrepreneurial Studies, Niigata, 950-0916 Japan
jma11202097zp@jigyo.ac.jp, sugimoto.hitoshi@jigyo.ac.jp, otsuka.akira@jigyo.ac.jp

Abstract. In China, "Green" is a priority of the national economic and social development
plan, and environmentally friendly manufacturing is important. In order to achieve environmentally friendly manufacturing, it is also necessary to optimize the allocation of energy resources in material procurement, production and distribution. In this paper, we willconsider the concept of environmentally friendly manufacturing.

Keywords: Green manufacture, environment-friendly

1. はじめに
持続可能な社会実現に向けて、環境への配慮が重要であり、製造業における環境の対応は喫緊の課題となっている。中国では、今後の国民経済・社会発展計画において、「グリーン（環境）」が重点事業として挙げられている。2021 年 3 月の全国人民代表
大会（全人代）で承認された「国民経済・社会発展第 14 次 5 ヵ年規画と 2035 年までの長期目標要綱」では、「グリーン発展の推進」が重要事項として今後 5 年間の活動内容として挙げられており、中国の発展において、イノベーションに加え、人と自然の調和が方策であることが示されている。[1], [2]中国において、イノベーションを進めるためには製造におけるグリーン化、「グリーン製造」が重要である。グリーン製造の実現のためには、製造機器の高度化に加え、65

JSST2021
資材調達、生産、物流におけるエネルギーなどの資源の投入の最適配置も必要となる。
本稿では、環境に配慮した製造である、グリーン製造の考え方について検討する。

2. 中国の「グリーン発展の推進」
中国第 14 次五ヵ年規画におけるグリーン発展の推進の概要を説明する。第 13 期全国人民代表大会（全人代）第 4 回会議は 2021 年 3 月に北京市で開幕され、第 14 次 5ヵ年（2021～2025 年）と 2035 年までの長期目標綱要案が発表された。第 14 次 5 ヵ
年計画におけるグリーン発展の推進の位置付けを図１に概要を示す。[3], [4]図 1 中国グリーン発展の推進の位置付け

第 14 次 5 カ年規画では、次の通り計画されている。第 11 篇 「グリーン発展の推進」は、人間と自然の協和共存を促進する。第 39 章「開発モードのグリーン化を加速する」は、生態の優先順位とグリーン開発を順守し、総合的な資源管理、科学的な配分、節約を行う。更に、リサイクルを促進し、質の高い経済発展と生態環境の高度の保護を促進する。
第 39 章の中の、第１節は「資源利用効率を総合的に向上」、第２節は、「資源循環システムを構築」、第３節は、「グリーン経済を積極的に発展させる」、第４節は、「グ
リーン開発政策体系の構築」である。グリーン製造関連の具体的な標準では、高エネ
ルギー、高排出なプロジェクトを無条件で実施することを抑制し、エコへの転換を推進する。また、省エネ、クリーン生産、クリーンエネルギー、カーボンニュートラル、

66

JSST2021
エアコンや暖房などの基礎的な設備をグリーン化すること、および電気自動車などのグリーンサービス産業を強化し、契約エネルギー管理、契約節水管理、環境汚染の管理を推進、グリーン発展の法律と政策保障を強化する。省エネ・環境保全と資源総合利用に有利な税収政策を実施する。また、グリーン化を進める企業への融資を優先する取り組みであるグリーン投資へも力を入れる。但し、自然資源の有償利用制度を健全化し、自然資源、汚水ごみ処理、水資源、エネルギーなどの分野では、各業界においては、自然資源の保護に努めるように価格を改善している。固定資産投資プロジェ
クトの省エネ審査、省エネ監察を行い、各業界のエネルギーの使用についてコントロールする管理制度の改革を推進する。

3. グリーン製造の考え
製造業のグリーン化においては。資材調達、生産、流通に加え、消費するエネルギ
ー、環境保護措置を考えることが必要である。まずは、必要な時間を必要な場所に投資する」というジャストインタイム(JIT)の考え方[5]を基本として、環境＆エネルギー（Environment&Energy）に配慮した生産を検討する。
図 2 に環境&エネルギージャストインタイム(E-JIT)の概要を示す[6]。

図2E-JIT 概要ある製品を製造する上では、生産資産（部品）の調達、および生産資源（電力、水など）が必要となる。そして、各生産工程で必要な部品、資源の配置、および生産工67

JSST2021
程で取られる環境保護措置を検討する必要がある。さらに、生産資産を生産現場へ運
ぶための資源も検討が必要となる。今後の製造においては、各生産工程と、生産資産、生産資源、環境保護措置、部品の物流の最適化が課題となる。
グリーン製造の評価においては、使用される生産資源、環境保護措置を最小化しつつ、製品製造を最大化する最適化の E-JIT 生産管理ツールが必要となる。今後は、EJIT 生産管理ツールの要件を検討する。

4. おわりに
本稿では、製造業における環境の対応として、中国におけるグリーン製造の考えとして E-JIT を検討した。今後、環境と製造の定量化のためのツールの検討を進め、中国グリーン製造の評価手法を構築する。

References[1] 日本貿易振興機構(JETRO), “全人代、次期 5 カ年規画と 2035 年までの長期目標綱要案を発表”, https://www.jetro.go.jp/biznews/2021/03/38a0764930a67b96.html(2021 年 6 月 21 日参照)

[2] 独立行政方針経済産業研究所, “始動する中国における第 14 次五ヵ年計画 ―「 質 の 高 い 発 展 」 を 目 指 し て ― ”, https://www.rieti.go.jp/users/chinatr/jp/210415kaikaku.html (2021 年 6 月 21 日参照)
[3] 中华人民共和国国家发展和改革委员会, ”中华人民共和国国民经济和社会发展第十四个五年规划和 2035 年远景目标纲要”
https://www.ndrc.gov.cn/xxgk/zcfb/ghwb/202103/t20210323_1270124.html (2021
年 6 月 21 日参照)
[4] THE STATE COUNCIL THE PEOPLE'S
lish.www.gov.cn/ (2021 年 6 月 21 日参照)

REPUBLICOF

CHINA,

http://eng-[5] 大野耐一, “トヨタ生産方式”, ダイヤモンド社, 1978 年 5 月[6] 富澤克行, ”三菱電機事業概要＆智能製造理念”, 瀋陽市智能製造解決方案説
明会, 2019 年 3 月

68

JSST 2021

NS3 を用いた無線通信規格が過密な Wi-Fi 環境に与える影響に関する研究
The study on effects of wireless standards ondense Wi-Fi environments with NS3
岩田 寛司 1*, 堀内

咲江 2

Hiroshi Iwata1*, Sakie Horiuchi2

1

岐阜工業高等専門学校 先端融合開発専攻
2

岐阜工業高等専門学校 電気情報工学科

Gifu National College of Technology, Department of Advanced Fusion Development1

Gifu National College of Technology, Department of Electrical and Information Engineering2*2021y10@stu.gifu-nct.ac.jp

Abstract. Currently, there are many wireless LANs. Therefore, some wireless LANs have no
choice but to use the same frequency channel, and it causes the degradation of throughput. In
this paper, we investigate the effects of wireless standards on dense wireless LANs environments. In particular, we focus on the throughput of wireless LANs based on IEEE802.11n,
IEEE802.11ac, and IEEE802.11ax. Further we analyze the throughput characteristics.
Keywords: Wireless LAN, Troughput, Capture effect, Interference, NS3

1. はじめに
現在の社会では，モバイル端末の利用増加に伴い，Wi-Fi も広く普及している．都
市部の駅やカフェなどの人の多く集まる場所では，モバイル Wi-Fi やテザリングを使用するユーザが増加し，Wi-Fi の密集地帯が発生している現状がある．そのような状況により，同じチャネルを使用する Wi-Fi が増加し，通信品質の劣化につながると懸念される．研究[1]では，過密な無線 LAN 環境に対して，干渉とキャプチャ効果を考慮したスループット計算シュミレータを開発した．これに加えて実験[2]では，同様な環境に対して実機を用いた実験により，各無線 LAN および全体のスループット特性について調査を行っている．これらの研究では，1 台のアクセスポイント（以下，AP）

69

JSST 2021に 1 台の端末が帰属したものを 1 つの無線 LAN とし，同一の周波数チャネルを使用する複数の無線 LAN が格子状に存在する環境でのスループット特性を調査している．
研究結果より，無線 LAN 同士の距離が 0.6m 程離れている場合，無線 LAN の台数が数台程度のときはキャプチャ効果によりスループットは向上するが，台数がさらに増加すると干渉の影響が強くなりスループットが劣化することが述べられている．しかしながら，これらの研究は無線通信規格として IEEE802.11g のみを対象としている．
近年では端末数の増加に対応した新たな無線通信規格が標準化されているため，これらの無線通信規格に対する干渉やスループット特性の究明が求められる．本研究では，
ネットワークシミュレータ NS3 を用いて無線通信規格が過密な Wi-Fi 環境に与える影響について調査を行う．限られた範囲内で同一の無線通信規格を使用する Wi-Fi を多数配置し，それによるスループット特性とその原因について考察していく．

2. 関連技術
近年，無線通信規格として，IEEE802.11n（以下，11n），IEEE802.11ac（以下，11ac），
IEEE802.11ax（以下，11ax）が標準化され，広く利用されている．11n は複数の送受信アンテナを搭載することで，2 つのデータストリームを同時に空間多重して送受信する MIMO 技術や，搬送波数の増加，隣接したチャネルを繋ぐチャネルボンディン
グにより，スループットの向上を実現する．11ac はチャネルボンディングや MIMOを拡張した MU-MIMO（Multi User Multi-Input Multi-Output）により，1 つのアクセス
ポイントから複数の端末に対して同時に送信することが可能となる．11ax では，複数のサブキャリアをユーザチャネルとして分割し，複数のユーザに割り当てることによって，多くのユーザの通信を同時に行うことができる変調方式である OFDMA や，1つのチャネルで 1 対多で通信を実現する MU-MIMO によって，多数の端末が同じチ
ャネルを使ってもスループットの減少を抑えられるようになっている．

3. 研究内容
本研究では，限られた範囲内で同一の無線通信規格を使用する Wi-Fi を多数配置し，通信を行った場合のスループット特性について調査する．本実験で使用するモデルを
図 1 に示す．1 台の AP と 1 台の端末のセットを 1 つの Wi-Fi とし，同一周波数のチャネルを使用する Wi-Fi を格子状に増加させた場合の総スループットの変化を調査する．この時，AP と端末の間には障害物はなく，Wi-Fi 同士の距離 d は，d=0.3，0.6[m]の値をとる．その他のパラメータは表 1 に示す．
各 Wi-Fi の使用する規格は統一とし，
11n，11ac，11ax のそれぞれの場合についてシミュレーションを行う．本研究ではネットワークシミュレータ NS3 を用いる．NS3 は離散イベントシミュレータであり，以下の式（1）で各 Wi-Fi のスループット Ti を導出する． Ti は送信されるパケット数
Pi からシミュレーション時間 S で割ることで導出できる．シミュレーション時間 S は，
全端末が未送信であった平均時間 S1，送信パケットの衝突が発生し，パケット送信が失敗した平均時間 S2，相手側でエラーが発生した平均時間 S3，パケット送信の衝突が起こらず正常に送信パケットが送れた平均時間 S4 の和である．本研究では，各 Wi-Fiのスループット Ti の総和をトータルスループットとする．
𝑃
𝑃𝑖
𝑇𝑖 ＝ 𝑆𝑖＝𝑆 +𝑆 +𝑆
…（1）
+𝑆
1

2

3

70

4

JSST 2021
表 1 パラメータ

図1

Traffic

飽和モード

伝送レート

Auto

通信プロトコル

UDP

トラフィックの方向

uplink使用チャネルCH1

Wi-Fi の配置

4. 数値例
本研究のシミュレーション結果を図 2，図 3 に示す．図 2 には Wi-Fi 間の距離
d=0.6[m]の場合，図 3 は d=0.3[m]の場合それぞれについて，配置する Wi-Fi の数を増加させた際のトータルスループットを示している．図 2 より，Wi-Fi 数が少数の場合
は，いずれの規格においても Wi-Fi が 1 台の時よりも高いスループットを示している．これはキャプチャ効果により，送信パケットの衝突が発生するにも関わらず，パケット送信が成功しスループットの向上に繋がったと考えられる．しかしながら，WiFi 数がさらに増加した場合，いずれの規格でもスループットが劣化しており，中でも11n のスループットの劣化が著しいことが分かる．理由として 11n では OFDM やMIMO などが用いられているため，干渉が大きくなった場合のスループットの劣化を抑えられなかったと考えられる．図 3 より，Wi-Fi 数が増加するといずれの規格においてもスループットが劣化している．とりわけ，11n と 11ac のスループットの劣化が顕著に現れており，一方で 11ax ではスループットの劣化が抑えられていることが分かる．これは，11ax の OFDMA によって多数の Wi-Fi が同じチャネルを使っても効率よく通信することが可能となり，スループットの劣化を抑えられたと考えられる．

図2

d=0.6[m]におけるスループットの変化

71

JSST 2021

図3

d=0.3[m]におけるスループットの変化

5. まとめと今後
本研究では NS3 を用いて限られた範囲内で同一の無線通信規格を使用する Wi-Fi
を多数配置し，それによるスループット特性を調査した．これにより，11n は干渉によるスループットの劣化が激しく，Wi-Fi 数が少数かつ密集していない環境でないと理想的なスループットを実現することができないことがわかった．11ac は Wi-Fi 間の距離が 0.6m 程度であれば，ある程度のスループットが担保されているが，それよりも近くに配置されるとスループットが急激に劣化した．一方で，11ax の Wi-Fi を近距離に配置した場合でもスループットの劣化を抑えられた．これによって，11ax は近距離に多数の Wi-Fi を配置してもスループットの劣化が抑えられると言える．今後は伝送レートの変化や再送率なども調査し，各規格にどのような特徴が現れるかをさらに検証していく予定である．

参考文献[1] 堀内咲江, 三好一徳, 村瀬勉, 橘拓至: 過密な無線 LAN におけるキャプチャ効
果を考慮したシミュレーション性能評価, 信学技報, vol. 114, no. 371, NS2014150 (2014), pp. 19-24[2] 堀内咲江, 村瀬勉, 橘拓至: 電波干渉とキャプチャ効果の影響に関する実験調
査, ITRC-NWGN(2017)72

JSST 2021

Computational method for interactions between
deformable objects and ﬂuid ﬂows using immersedboundary method and mass spring model
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

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.14.78

73

JSST2021

Simulation of elastic waves in micropolar bimaterials
using 2-D M-EFIT
Yusuke Suzuki1* and Takahiro Saitoh2
Department of Civil and Environmental Engineering, Gunma University, Japan

1, 2

*

t160c050@gunma-u.ac.jp

Abstract. In this research, a micropolar elastic wave simulation tool, which is called M-EFIT(Micropolar-Elastodynamic Finite Integration Technique), is developed. As numerical examples, elastic wave propagation and scattering in a micropolar bimaterial are demonstrated by
the proposed method.
Keywords: Elastic waves, Micropolar-elastodynamic finite integration technique (M-EFIT), Finite
difference time-domain method (FDTD), Micropolar-elastodynamics

1. Introduction
The concretes, bedrocks, and bones are known as microscopic materials. In general, the microscopic property is not considered in the classical elastodynamic theory. In order to integrate themicroscopic property into the classical elastodynamic theory, Eringen proposed the micropolar
elastodynamic theory. However, elastic wave simulation tools based on the micropolar elastodynamic theory have hardly been developed. Therefore, in this research, a micropolar elastic
wave simulation tool, which is called M-EFIT (Micropolar Elastodynamic Finite Integration
Technique), is developed. The EFIT is a kind of the FDTD (finite difference time-domain
method), and has been developed by Nakahata et al. [1]. This EFIT is extended to 2-D micropolar-elastodynamic problems in this research. Some numerical examples for wave propagationand scattering in a 2-D micropolar bimaterial are demonstrated by the proposed method.

2. M-EFIT formulation
The EFIT is a grid-based numerical simulation method based on the FDTD ,and can easily treat
the boundary conditions on the interface between different materials. Considering the 𝑥1 − 𝑥3
plane as a 2-D inplane problem, the equation of motions for the stress 𝜎𝛽𝛼 and the couple
stress 𝑚𝑖𝑗 , respectively, can be rewritten as follows:𝜌𝑢̈ 𝛼 =

𝜕𝜎𝛽𝛼
𝜕𝑥𝛽

𝜕𝑚
+ 𝜌𝑏𝛼 , 𝜌𝑗𝜙̈2 = 𝛼2 + 𝑒2𝛼𝛽 𝜎𝛼𝛽 + 𝑙2
𝜕𝑥𝛼

(1)

where 𝜌 is the density of a micropolar elastic material, 𝑢𝛼 is the displacement, 𝑗 is77

JSST2021
the micro moment of inertia, 𝜙2 is the micro rotation, 𝑏𝛼 is the body force, and 𝑙2 isthe body couple force. 𝑒2𝛼𝛽 is the permutation symbol and ( ̇ ) show the partial derivative with respect to time 𝑡. In addition, the constitutive equations for the stress
𝜎𝛼𝛽 and the couple stress 𝑚𝑖𝑗 can be obtained as follows:𝜎̇𝛼𝛽 = 𝜆

𝜕𝑢̇ 𝛽
𝜕𝑢̇ 𝛾
𝜕𝑢̇ 𝛼
𝛿𝛼𝛽 + (𝜇 + 𝜅)+𝜇
− (𝜇 + 𝜅)𝑒2𝛼𝛽 𝜙̇2 − 𝜇𝑒2𝛽𝛼 𝜙̇2
𝜕𝑥𝛾
𝜕𝑥𝛼
𝜕𝑥𝛽
𝑚̇𝛼2 = 4𝐺𝑁 2 𝑙 2

𝜕𝜙̇2
𝜕𝑥𝛼

(2)
(3)

where 𝛿𝛼𝛽 is the Kronecker delta, and the parameters, 𝜆, 𝜇, 𝜅, 𝐺, 𝑁 and 𝑙 are the micropolar elastic constants [2]. Equations (1)-(3) are discretized by using the EFIT
scheme for time and space [1]. The detail is skipped for the page limitation.

3. Numerical example
Numerical examples can be seen in Fig.1. In this simulation, the grid size 𝛥𝑑 and time step
size Δ𝑡 are given by 𝛥𝑑 = 1.0 × 10−2 and 𝛥𝑡 = 1.0 × 10−4 , respectively. The upper andlower materials are a micropolar and a conventional elastic material, respectively. The micropolar elastic constants are given by 𝜇 = 1.0, 𝜆 = 1.5, 𝜅 = 1.0 and 𝜌 = 4.5. On the other hand,
the lower conventional isotropic material parameters are 𝜇 = 1.0, 𝜆 = 1.5, 𝜅 = 0.0 and 𝜌 =4.5. A sinusoidal particle velocity is given at the top center of the upper material. As shown in
Fig.1(a) and (b), primary and secondary waves (P and S waves) exist in the both upper andlower materials, and reflected by the interface drawn by the white line in Fig.1. However, themicro rotational wave (M wave) can exist in the micropolar material only.

4. Conclusion
In this research, M-EFIT was developed and some numerical examples are demonstrated. In
the near future, this work will be extended to 3-D micropolar elastodynamic problems.

Fig.1: Elastic wave propagation in a micropolar bimaterial at (a) 80Δ𝑡 for P and S waves, (b)190Δ𝑡 for P and S waves (c) 80Δ𝑡 for M wave and (d) 190Δ𝑡 for M wave.

References[1]

K. Nakahata, S. Hirose, F. Schubert, B. K𝒐̈ hler: Image based EFIT simulation for nondestractive ultrasonic testing of austenitic steel, Journal of solid mechanics and materials engineering, Vol.3(12), (2009), 1256-1262.

[2]

T. Fukui, Y. Okui: Boundary element analysis of scattering problems in two dimensional
micropolar elasticity, Elsevier, 283-292.

78

JSST2021

Molecular Dynamics Study of Water Dynamics
in Zwitterionic Polymer Brush-Water Interface
Yuya Fujinaga1*, Susumu Fujiwara2, Tomoko Mizuguchi2
Graduate school of science and technology, Kyoto Institute of Technology
Faculty of Materials Science and Engineering, Kyoto Institute of Technology*m0672023@edu.kit.ac.jp
12

Abstract. Polymer brushes have a protein antifouling effect to a biomaterial in vivo. The dynamics of water in polymer brush-water interface is studied by the molecular dynamics simulations to clarify the molecular mechanism of the protein antifouling effect. In addition, we also
clarify the relationship between the structural characteristics of polymer brushes and the dynamics of water.
Keywords: Polymer brush, Molecular dynamics simulation, Dynamics of water, 2-methacryloyloxyethyl phosphorylcholine, Hydrogen bond, Tacticity

1.

Introduction

The surface properties of a material can be modified by polymer brushes which are obtainedby growing polymers on the surface of the material at high density. Polymer brushes give a
protein antifouling effect to the material[1]. The protein antifouling effect is very useful for
preventing protein adsorption which triggers the biological defense system when the material
is used in vivo, and therefore the application of polymer brushes to the medical field has attracted attention in recent years. As grafted polymer chains for polymer brushes, we selected
twitterionic 2-methacryloyloxyethyl phosphorylcholine (MPC), a biomaterial used in various
medical devices. The dynamics of water in polymer brush-water interface is important for understanding the protein antifouling effect. Clarification of the dynamics of water enables the
production of MPC polymer brushes with higher function, and is expected to lead further use
in the medical field. In the field of polymer brushes research, the type of a polymer and thecontrol of graft density have received much attention, but structural characteristics such as thetacticity of a polymer chain and graft position of polymer brushes have received less attention.
In this study, we focused on the relationship between structural characteristics and the dynamics of water. To this end, we performed molecular dynamics simulations of the systems composed of water and MPC polymer brushes with various structural characteristics.

2.

Simulation model and method

2.1 Moleculer model

79

JSST2021
We used an all-atom MPC polymer brush model in which
we grafted polymer brushes on the x-y plane (8 nm×8 nm)of the MD box. Water molecules were placed on the layerof the polymer brushes. As an initial configuration, weused the model after performing energy minimization(Fig.1). We prepared MPC polymer brush by using Winmostar[2]. As for force fields, we adopted TIP3P for waterand general AMBER force fields (GAFF) for MPC polymer brush[3]. As a point charge for each atom, we adopted
the Restrained Electrostatic Potential (RESP) charge calculated by Winmostar.
2.2 Simulation method
We performed MD simulations by Gromacs[4] with themodel under the NPT (constant number of atoms, constant pressure, and constant temperature) emsemble and2D periodic boundary conditions in x and y directions.
Since this study assumes the condition in vivo, we set T =310 K and P = 1 atm.

3.

Figure 1 Initial configuration of oursimulation model. Polymer brushes areplaced below water molecules.

(a)Simulation Results

We show snapshots of polymer brush-water interface at 0
ps and 10 ps in Fig. 2. We can see from this figure thatwater molecules penetrated into the layer of the polymer
brushes at 10 ps while the water layer and the polymerbrush layer are separated in the initial configuration at 0
ps. From our simulation data, we calculate the number
and the lifetime of hydrogen bonds, the diffusion coefficient of water, and so on, to understand the dynamics ofwater in polymer brush-water interface.

(b)

Figure 2 Snapshots of polymer
brush-water interface (a) at 0 ps and(b) at 10 ps.

References[1] Y. L. Liu, D. Zhang, B. P. Ren, X. Gong, L. J. Xu, F. A. Zhang, Y. Chang, Y. He, and J.
Zheng: Molecular simulations and understanding of antifouling zwitterionic polymer
brushes, J. Mater. Chem. B, 8:17(2020), 3814-3828.
[2] https://winmostar.com/jp/.
[3] topolbuild1_2_1.tgz in
http://www.gromacs.org/Downloads/User_contributions/Other_software.
[4] https://www.gromacs.org/.

80

JSST 2021

高周波電磁界ソルバ：
ADVENTURE_FullWave の高速化開発
Development for High-Speed of High-Frequency
Electromagnetic Fields Solver:
ADVENTURE_FullWave
大中 健登 1*
Kento Ohnaka1*
武居 周 1
Amane Takei1

1

宮崎大学 工学部 電気システム工学科

1

Department of Electorical and Systems Engineering, Faculty of Engineering,Miyazaki University*

hl16004@student.miyazaki-u.ac.jp

Abstract. Thigh-frequency electromagnetic fields solver: ADVENTURE_FullWave is a finite element analysis software in the parallel computational mechanics system: ADVENTURE System. The ADVENTURE_FullWave can perform detailed, high-speed, andhigh-efficiency finite element analyses of large-scale electromagnetic field problems using
the hierarchical domain decomposition method (HDDM) and the corresponding parallel distributed processing environment. The result of performance evaluations of the ADVENTURE_Fullwave will be shown in the conference.
Keywords: Finite element method, Electromagnetic fields, Large-scale analysis, Hierarchicaldomain decomposition method

1. 背景
近年の計算機の性能向上と数値計算技術の進歩により，計算機を用いた，高周波帯域の電磁界解析が普及した．これにより，アンテナ等の開発・設計において，実際に機器を試作し，測定を行わなくとも，電磁界の振る舞いを予測することが出来るようになり，開発期間の短縮や開発コストの削減に非常に有用なツールとして広く利用されるようになった．オープンソースソフトウエア ADVENTURE_FullWave[1]もその１つであり，解析モデルの作成・メッシュ生成・境界条件設定などで ADVENTURE

81

JSST 2021
System の各種モジュール群を基盤に用いた解析システムの構築とともに運用される．
ADVENTURE_FullWave は大規模解析に向けた並列計算手法として領域分割法を適用しており，また領域分割法に基づく並列化手法として COCG 法などの反復法を用いた反復型領域分割法を，さらに階層型領域分割法（HDDM）の適用により複数の計算機による並列計算・負荷分散を実現している ．これに加えて基盤とする
ADVENTURE System による特徴として，数千万～数億に及ぶ自由度メッシュによる詳細な解析が行えること，数万に及ぶプロセッサから成る超並列計算環境においても高い並列効率を得られること，PC クラスタから超並列計算機まで多様な並列分散環境に対して優れた移植性をもつことなどが挙げられる．
このように本システムには多様な並列分散環境に対応しながら高精細モデルを用いた大規模解析が行える利点がある一方，計算時間がかかるという高周波解析特有の問題があるため，特に並列計算におけるハイブリッド並列化の適用や線形代数ソルバの変更などの，並列計算機環境上でのコードチューニングにより処理速度のさらなる高速化や反復法の収束性改善が求められる．
そこで，本研究ではまずハイブリッド並列化をはじめとするコードチューニングによる高速化を検討し，今後の開発指針を得る．

2. 基礎式
2.1. ベクトル波動方程式
解析領域 Ω を考え，その境界を∂Ω とする．変位電流を含む Maxwell 方程式より導かれる電界 E[V/m]を未知関数としたベクトル波動方程式(1)を考える．
rot {(1 ⁄ 𝜇) rot 𝐄} − ω^2 ε𝐄 = jω𝐉
𝐄 ×𝐧 = 0

in Ω

(1a)
(1b)

on ∂Ω

(1c)

𝐉 = 𝜎𝐄̂

電界は単一角周波数 ω[rad/s]をもち，μ は透磁率[H/m]，ε は複素誘電率[F/m]，σは導電率[S/m]，J は電流密度[A/m2]を表す．ｊは虚数単位（ j = √−1 ），n は境界における外向き法線ベクトルである．複素誘電率 ε は ε = −𝜎/jω で与えられる．解くべき方程式は式(1a)であり，電界 E を未知数とした E 法として定式化される．
本定式化において単一角周波数 ω，透磁率 μ，複素誘電率 ε を既知量とし，放射源として電流密度 J を直接，あるいは式(1c)より既知電界 𝐄̂ を与える．また，得られた電界 E より Maxwell 方程式の１つである式(2)を用いて磁界 H[A/m]を計算する．
rot 𝐄 − jω𝜇𝐇 = 0(2)2.2. 有限要素定式化
式(1a)より弱形式を導き，領域 Ω の有限要素分割を考える．電界 E を Nedelec 四面体一次要素（辺要素）で近似し，電流密度 J を通常の四面体一次要素で近似する．
𝐄h , 𝐉h をそれぞれ E , J の有限要素近似とする．𝐄h∗ は任意の試験関数とする．

82

JSST 2021
2

∭ (1⁄𝜇) rot 𝐄h ∙ rot 𝐄h∗ dv − ω ∭ ε𝐄h ∙ 𝐄h∗ dv = jω ∭ 𝐉h ∙ 𝐄h∗ dv
Ω

Ω

(3)

Ω

式(3)を解くべき有限要素方程式とする．

3. 階層型領域分割法（HDDM）
並列有限要素法において最も有効な解法の１つであるとされる領域分割法では，解析領域全体をいくつかの部分領域に分割し，各部分領域に関する有限要素計算と，各部分領域境界上のつり合いをとる反復計算とを交互に行い，部分領域境界のつり合いが取れた時点で領域全体の解とする．解析領域の自由度は部分領域内部の自由度と境界上の自由度とに分けられ，静的縮約により内部の自由度を消去し境界上の自由度のみからなる方程式（インターフェース問題）を導出する．インターフェース問題を解く方法として共役勾配法（CG 法）などの反復法を用いる反復型領域分割法を用いる
が，CG ループ内各ステップでの部分領域内部の自由度を求める箇所において計算負荷が特に大きくなることがわかっている．この増大負荷を超並列計算機上で調整しながら効率的に解くための手法として階層型領域分割法を導入する．
図 1 および 2 に階層型領域分割法における領域分割とデータ分配の例を表す模式図を示す．階層型領域分割法では，まず全体領域を複数のパート（part）に分割し，さらに各 part をいくつかの部分領域（subdomain）に分割することで２段階の階層をもつ領域分割を行う．part 単位のデータは同個数の Parent と呼ばれるプロセッサに分配され，メモリ上に記憶される．これにより Parent の数を増やすほどより大規模な計算に対応することが可能となる．計算負荷が大きく増大する部分領域内部の自由度を求める箇所において，subdomain を適切な大きさ（ADVENTURE_FullWave
では１subdomain あたり要素数にして 100~200 が適当）に調節して分割していくことで，効率的な計算処理を実現できる．

図１．階層型領域分割法による領域分割

83

JSST 2021

図２．複数プロセッサへのデータ分配

4. 計算機環境および性能評価
性能評価に用いる計算機は，筆者らの研究室で所有する PC クラスタ Tomahawk である．本計算機は，Intel Core i7-2600K のマルチコア CPU および 32GB のメモリが搭載されている PC クラスタであり，25 台(100 core)を用いた．用いたコンパイラは gcc である．また，並列化ライブラリとして MPI(Message Passing Interface)を用いる．また，ダイポール・アンテナモデル[3]を用いた計測結果を表 1 に示す．本
結果に関する詳細な考察は，講演にて述べる．
表 1.計測結果
解析周波数

1.0 GHz

要素数

15,712,684

辺自由度数

5,599,463

計算時間[s]1,998

反復回数[回]

6,578

参考文献
[1] 武居 周, 室谷 浩平, 吉村 忍, 金山 寛: 数値人体モデルを用いたマイクロ波帯域の有限要素電磁界解析, 日本シミュレーション学会論文, Vol.4:No.3 (2021),
pp.81-95[2] 武居 周, 杉本 振一郎, 荻野 正雄, 吉村 忍, 金山 寛: 階層型領域分割法において部分領域に直接法を適用した高周波電磁場の大規模解析, 電気学会論文誌 A（基礎・材料・共通部門誌）, Vol.130:No.3 (2010), pp.239-246[3] 東 勇, 山田 知典, 武居 周: 反復型領域分割法に基づくマイクロ波シミュレー
ション, 日本シミュレーション学会論文誌, Vol.10:No.1 (2018), pp.25-33

84

JSST 2021

波動音響解析ソルバ：
ADVENTURE_Sound の開発と応用

Development and application for wave-sound field
Solver: ADVENTURE_Sound
吉留 涼 1*
Ryo Yoshidome1*
武居 周 1
Amane Takei1

11

宮崎大学 工学部 電気システム工学科

Department of Electorical and Systems Engineering, Faculty of Engineering,Miyazaki University*

hl18048@student.miyazaki-u.ac.jp

Abstract. The wave-sound fields analysis solver: ADVENTURE_Sound is a finite element
analysis software in the parallel computational mechanics system: ADVENTURE System. The
ADVENTURE_Sound can perform detailed, high-speed, and high-efficiency finite element
analyses of large-scale wave-sound field problems using the hierarchical domain decomposition method (HDDM) and the corresponding parallel distributed processing environment. The
result of performance evaluations of the ADVENTURE_Sound will be shown in the conference.
Keywords: Finite element method, Wave-sound fields, Large-scale analysis, Hierarchical domain
decomposition method

1. 背景
近年の計算機の性能向上と数値計算技術の進歩により，計算機を用いた，波動音響解析が普及した．これにより，コンサートホールや楽器の設計において，実際に実物やモックアップを試作し，測定を行わなくとも，音の振る舞いを予測することが出来るようになり，コンサートホールの建設や楽器の制作期間の短縮やコストの削減に非常に有用なツールとして広く利用されるようになった．オープンソースソフトウエア
ADVENTURE_Sound[1]もその１つである．ADVENTURE_Sound は, これまで小規模に留まっていた解析規模を一気にライブハウス・コンサートホールのスケールの計算を可能とした並列計算コードである．並列化手法として，領域分割法を適用しており，また領域分割法に基づく並列化手法として COCG 法などの反復法を用いた反85

JSST 2021
復型領域分割法を，さらに階層型領域分割法（HDDM）の適用により複数の計算機による並列計算・負荷分散を実現している．加えて基盤とする ADVENTURE System
による特徴として，数千万～数億に及ぶ自由度メッシュによる詳細な解析が行えること，数万に及ぶプロセッサから成る超並列計算環境においても高い並列効率を得られること，PC クラスタから超並列計算機まで多様な並列分散環境に対して優れた移植性をもつことなどが挙げられる．
本研究では，ADVENTURE_Sound を数千万要素を超える高精細モデルによる解析が必要となるバイオリンの設計に適用するために，必要となる並列計算機能の強化として，まずハイブリッド並列化をはじめとするコードチューニングによる高速化を検討し，今後の開発指針を得る．

2. 基礎式
2.1. Helmhortz 方程式
指定境界内部の解析領域Ωにおいて，音源をもつ音場内の速度ポテンシャルに関する波動方程式は次式で表される[2]．
𝛻!𝜙 −

1 𝜕!
𝜙=𝑞
𝑐 ! 𝜕𝑡 !

(1)

ただし，𝜙は速度ポテンシャル，cは音速，qは音源の分布関数である．
ここで，定常状態を考えるため，𝜙 = 𝛷𝑒𝑥𝑝 (−𝑗𝜔𝑡)とすると(1)式は
𝛻!𝛷 +

𝜔!
𝛷=𝑞
𝑐!

(2)

となり，解くべきHelmholtz方程式を得る．ここに，𝑗は虚数，𝜔は角周波数である．
また，音圧𝑝は未知数である速度ポテンシャル𝛷を求めたのちに次式で求められる．
𝑝 = −𝑗𝜔𝜌𝛷(3)
3

ここに，𝜌は媒質の密度[kg/m ]である．

2.2. 有限要素定式化
式(2)にガラーキン法を適用することで弱形式化し，有限要素法により近似すると次式が得られる[3]．
! ∇𝛷! ∙ ∇𝛷! ∗ 𝑑𝛺# −
$!

𝑗𝜔𝜌
2
+ 𝛷 𝛷 ∗ 𝑑𝛤 − 𝑘 ! 𝛷! 𝛷! ∗ 𝑑𝛺# = 0.
𝑍% &! ! ! #𝛺(4)(4)式を解くべき有限要素方程式とする．𝛷は未知数となる速度ポテンシャル， 𝑘と𝜔は波数ならびに角周波数, 𝜌は密度, そして𝑍! 固有音響インピーダンスである.

86

JSST 2021

3. 計算機環境および性能評価
性能評価に用いる計算機は，筆者らの研究室で所有する PC クラスタ Tomahawk である．本計算機は，Intel Core i7-2600K のマルチコア CPU および 32GB のメモリが搭載されている PC クラスタであり，25 台(100 core)を用いた．用いたコンパイラは gcc である．また，並列化ライブラリとして MPI(Message Passing Interface)を用いる．また，バイオリン簡易モデルを用いた計測結果を表 1 に示す．本結果に関する詳細な考察は，講演にて述べる．
表 1.計測結果
解析周波数[Hz]

500Hz

要素数

2,956,804

自由度数

4,005,344

計算時間[s]23.59s

反復回数[回]

14回

参考文献[1] ADVENTURE Project Homepage [Online]: https://adventure.sys.t.u-tokyo.ac.jp/ .
[2] N. Sakamoto, T. Otsuru, R. Tomiku, S. Yamauchi, “Reproducibility of sound absorption and surface impedance of materia
ls measured in a reverberation
room using ensemble averaging technique with a pressure-velocity sensor andimproved calibration,” Applied Acoustics, vol.142, no.15, pp.87-94, 2018.
[3] 武居 周, 工藤彰洋, “並列有限要素法に基づく 1 億自由度超の波動音響解析,” 日
本シミュレーション学会論文誌，vol.12, no.2, pp.76-84, 2020.

87

JSST2021

Molecular dynamics simulation and numeric calculation of a damaged polyethylene assuming tritium substitution and decay
Ryuta Kawanami1,2*, Susumu Fujiwara3, Hiroaki Nakamura4,5, Kazumi Omata6,7
Graduate School of Science and Technology, Kyoto Institute of Technology
Quantum Beam Science Research Directorate and Institute for Quantum Life Science,
National institute for Quantum and Radiological Science and Technology
3
Faculty of Materials Science and Engineering, Kyoto Institute of Technology
4
Department of Helical Plasma Research, National Institute for Fusion Science5
Department of Energy Engineering and Science, Nagoya University
6Clinical Research Center, National Center for Global Health and Medicine
7Departments of Hematology, Rheumatology, and Infectious Diseases, Kumamoto University
1

2*

m0672007@edu.kit.ac.jp

Abstract. It is important to simulate how the structure of a molecule, such as DNA, responds
when subjected to tritium substitution and decay. This is because structural changes in DNA
can inhibit enzymatic repair. We have observed the response of single-stranded polyethylene,
a simpler polymer, by molecular dynamics simulations and attempted to make theoretical predictions based on non-equilibrium statistical mechanics. As a result, we found that we could
predict the response to some extent for low damage rates.
Keywords: Molecular dynamics simulation, Linear response theory, Theoretical calculation, Tritium, Polymer, Decay effect

1. Introduction
The accident at the Fukushima Daiichi Nuclear Power Plant in 2011 resulted in the generationof a large amount of radiation-contaminated water. Although the radionuclide removal system(ALPS) can remove radioactive materials other than tritium, the treatment of the resulting tritiated water has become a problem. The main target of tritium is DNA, which is responsible forthe production of proteins. Although the effects of tritium on macromolecules have been studied [1,2,3], the biological effects of tritium have not yet been elucidated. One possibility is thatthe structural changes in the damaged DNA may inhibit enzymatic repair. There is also thepossibility that strand breaks will be induced. Therefore, it is important to discuss the structuralchanges of tritium-damaged DNA. However, many interactions such as base stucking, hydrogen bonding between interstrand bases, and shielding of the negative charge of phosphate88

JSST2021
groups by cations contribute complexly to the stability of DNA. Therefore, in this study, we
will investigate the structural changes in polyethylene, a simpler polymer macromolecule, instead of those in DNA, a complex macromolecule.
Tritium beta decay occurs with a short half-life of 12.3 years, producing helium-3. The rangeof the beta-electrons in water is at most 6m, so external exposure is not a problem. What canbe a problem is when tritiated water is absorbed into the human body; either the tritiated water
near the DNA emits beta electrons, or the hydrogen in the DNA in contact with the tritiatedwater is replaced by tritium and it decay. We focused on the latter, the decay of tritium in DNA
given by tritiated water and the change in its structure. This is a unique effect of tritium, aradioactive isotope of hydrogen.
The structural change after tritium damage can be studied by molecular dynamics (MD) simulations. Another approach is to predict the response function using theories in nonequilibrium
statistical mechanics. It has been suggested that this can be used to calculate the response ofdamaged molecules from MD simulations of equilibrium (undamaged) system [4]. In the second half of this paper, we will compare the response functions calsulated by the theory withactual MD simulations of damaged molecules and discuss their validity.

2. Simulation
In this study, a model of unbranched, single-chain polyethylene was used. Each carbon andthe hydrogen bonded to it is coarse-grained, with a particle number of 3000. Damage due totritium substitution and decay was represented by changing the type of coarse-grained particle
to one with one hydrogen removed (e.g. CH2 → CH, CH3 → CH2). Assuming that the damage
caused by tritium is completely random, the ratio for the number of particles is fH.
Try to evaluate the structural change of polyethylene with the global orientation order parameter 𝑃2 as follows.
3 cos2 𝛷 − 1
𝑃2 = ⟨⟩
3all−bond

When polyethylene is at 250K in vacuum (Figure 1), the global orientation order parameter
P2 of polyethylene with 𝑓H = 0 is 0.431. On the other hand, P2 for 𝑓H = 0.1 is 0.404 andthere is no significant difference. We decided to investigate the potential energy of the system.

Figure 1. Top and side viewsof polyethylene; 𝑓H = 0 (left),
𝑓H = 0.1 (right). The large
spheres indicate the damagedsites.

Figure 2. Potential energy relaxation of damaged
polyethylene. They are approaching the value of𝑓H = 0 with oscillation.

89

JSST2021
Figure 2 shows the relaxation of the potential energy 𝐸 = 𝐸bond + 𝐸angle + 𝐸torsion + 𝐸LJof the system at 𝑓H = 0.1 and 𝑓H = 0.01 at 250K. The value settled at equilibrium is thesame as that for 𝑓H = 0, which means that the polyethylene transitions to a slightly different
structure after the decay. In fig. 1, the potential energy oscillates with a buzzing motion, whichis thought to be the oscillatory particles with different frequencies or masses.

3. Linear response theory
Linear response theory is used to predict the response of damaged polyethylene. In this theory,
the response function of a system subjected to a weak perturbation is obtained as a functionsof the physical quantity in equilibrium. Suppose that in the past when 𝑡 < 0, the system follows the Hamiltonian 𝐻 0 . When a perturbation 𝐻 ′ (𝑡) is added at the moment 𝑡 = 0 , thephysical quantity 𝐵 of the system evolves according to 𝐻(𝑡) = 𝐻 0 + 𝐻 ′ (𝑡). Suppose that the
perturbation term 𝐻 ′ (𝑡) of the Hamiltonian can be devided into 𝐹(𝑡) and 𝐴(𝛤) , where𝐹(𝑡) is a time-dependent external field term and the conjugate 𝐴(𝛤) is a quantity that depends on the phase spase. In this study, we consider tritium damage as a perturbation, which is
given by 𝑓H simultaneously at 𝑡 = 0. The physical quantity we want to predict is the potential
energy. That is, let 𝐹(𝑡) be a unit step function with value at 𝑡 > 0 and 𝐴(𝛤) be the potential energy of the system itself. Linear response theory gives the response function 𝜙(𝑡), which
is a function only of 𝐵 in equilibrium and the time derivative 𝐴̇ of the potential energy. Specifically, it is expressed as a time correlation function of 𝐵 and 𝐴̇ at equilibrium.
(1)
𝜙𝐵𝐴 (𝑡) = 𝛽⟨𝐵(𝑡)𝐴̇(0)⟩𝑡⟨Δ𝐵⟩ = ∫ 𝑑𝑡 ′ 𝜙𝐵𝐴 (𝑡 − 𝑡 ′ )𝐹(𝑡 ′ )

(2)

−∞

4. Discussion
We first analytically derived the 𝐴̇ corresponding to 𝐸bond , 𝐸angle and 𝐸torsion , respectively. Then, we substituted the result of 𝑓H = 0 into that equation and calculated the response
function numerically for 𝑓H = 0.1. The ensemble average in Equation (1) is for 2000 samples.
In the integration of the response function in Equation (2), the calculated response function isshifted so that the values do not diverge.
The result show that for 𝑓H = 0.1 , thepredictions differ by a factor of 10 (Figures
1 and 2). A possible reason for this is thatfor a large fraction of damage, 𝑓H = 0.1,
the nonlinear term is large, i.e., the damagesite has a strong effect of changing thestructure in a concerted manner.
In the future, we will verify whether the
results are consistent with the actual simulation values for 𝑓H = 0.01 and even Figure 3. Relaxation of polyethylene potensmaller fractions. We will also calculate the tial energy calculated based on linear reglobal orientation order parameter 𝑃2 and sponse theory.
the other quantities that represent structuralchanges as physical quantities 𝐵 to be predicted.

90

JSST2021

References[1] S. Fujiwara, H. Nakamura, H. Li, H. Miyanishi, T. Mizuguchi, T. Yasunaga, T. Otsuka,
Y. Hatano, S. Saito: Computational strategy for studying structural change of tritiumsubstituted macromolecules by a beta decay to helium-3, J. Adv. Simulat. Sci. Eng., 6:1(2019) 94-99.
[2] H. Li, S. Fujiwara, H. Nakamura, T. Mizuguchi, T. Yasunaga, T. Otsuka, T. Kenmotsu,
Y. Hatano, S. Saito: Structural changes in tritium-substituted polymeric materials by abeta decays: A molecular dynamics study, Plasma Fusion Res., 14 (2019) 3401106.
[3] H. Li, S. Fujiwara, H. Nakamura, T. Mizuguchi, A. Nakata, T. Miyazaki, S. Saito: Structural change of damaged polyethylene by beta-decay of substituted tritium using reactive
force field, Jpn. J. Appl. Phys., 60 (2021) SAAB06.
[4] R. Kubo: Statistical-Mechanical Theory of Irreversible Process. I. General Theory andSimple Applications to Magnetic and Conduction Problems, J. Phys. Soc. Jpn., 12:6
(1957), 570-586

91

JSST 2021

void 欠陥を有する磁性体の磁化過程解析
Magnetization process analysis of magnetic material with void defect五十嵐靖幸*, 山口克彦
Yasuyuki Igarashi*, Katsuhiko Yamaguchi

福島大学 共生システム理工学研究科
Graduate School of Symbiotic Systems Science and Technology, Fukushima University*

s2170005@ipc.fukushima-u.ac.jp

Abstract. Assuming iron clusters with void defects, we analyzed how the defects affect themagnetization process using the LLG equation. When the cluster had void defects, changes ofenergy were seen in the exchange energy and magnetic dipole energy as the domain wall
passed through the defect
Keywords: LLG equation, micro magnetic simulation, domain wall displacement, void defect1. はじめにバルクハウゼンノイズ(BHN)は, 鉄などの磁性体が磁化される過程で, 磁壁が不連続に移動することで発生する [1]. 磁壁の不連続な移動は磁性体内の欠陥や析出物などに, 磁壁移動がトラップされることによって引き起こされる[1]. すなわち BHN は
磁性体内の欠陥に敏感であることから, 鉄鋼材料の劣化具合を測定する非破壊検査に用いられている. 磁壁が磁性体内の void 欠陥を通過する際に, 磁壁は欠陥にトラップされて磁気双極子エネルギーが減少する説が唱えられている[2]が, 欠陥が磁壁に与える影響を微視的に解析した例は少ない. そこで微視的な磁化解析に用いられる
Landau-Lifshitz-Gilbert(LLG)方程式を使用して, void 欠陥を含む磁性体の磁化過程解析をおこなった.

2. 手法
まず磁化の時間変化を計算する際に用いられる LLG 方程式について(1)式に示す.
𝑑𝑴⁄𝑑𝑡 = −𝛾𝑴 × 𝑯𝒆𝒇𝒇 + (𝛼⁄𝑀𝑠 )𝑴 × 𝑑𝑴⁄𝑑𝑡 (1)(1)式中の𝑴は磁化ベクトル, γはジャイロ磁気定数(𝛾 = 1.76 × 107 [1/𝑂𝑒 ⋅ 𝑠]), 𝑯𝒆𝒇𝒇 は有効磁場, αはギルバート減衰定数(0 ≤ 𝛼 ≤ 1), 𝑀𝑠 は自発磁化である.

92

JSST 2021
𝑯𝒆𝒇𝒇 を算出する式を(2)式に示す.
𝑯𝒆𝒇𝒇 = − ∂𝐸𝑡𝑜𝑡 ⁄∂(𝑀𝑠 𝛼)(2)
また, 𝐸𝑡𝑜𝑡 について(3)式に示す.
𝐸𝑡𝑜𝑡 = 𝐸𝑒𝑥𝑡 + 𝐸𝑎𝑛𝑖 + 𝐸𝑒𝑥 + 𝐸𝑑𝑝(3)
𝐸𝑡𝑜𝑡 は以下で述べるエネルギーの総和である. ゼーマンエネルギー𝐸𝑒𝑥𝑡 は外部から磁性体中に印加した磁場によるエネルギー, 異方性エネルギー𝐸𝑎𝑛𝑖 は磁性体中の磁気モーメントを磁化しやすい方向(容易軸)から磁化しにくい方向(困難軸)にねじ向けるのに必要なエネルギー, 交換エネルギー𝐸𝑒𝑥 は磁性体中で隣接する磁気モーメントに働
く交換相互作用によるエネルギーで, 磁気双極子エネルギー𝐸𝑑𝑝 は磁性体中の磁気モーメント同士に働く磁気双極子相互作用によるエネルギーである.これらのエネルギ
ーは磁気モーメントの方向より算出させる.
計算には, M.R.Scheinfein らによって作成された LLG Micromagnetic Simulator(以下,LLG シミュレータ)という LLG 方程式を用いて解析を行えるシミュレーションソフト
を用いた. LLG シミュレータでは, 計算を想定する磁性体クラスタの大きさ, 想定する材料の物性パラメータ, 境界条件, 計算開始時の磁性体の磁化(初期磁化), 印加する外部磁場, 欠陥の設定をすることができる. 物性パラメータには材料ごとに異なるものとして, 𝑀𝑠 , 交換スティフネス定数𝐴, 異方性定数𝐾がある. 𝐴は交換エネルギーの計算, 𝐾は異方性エネルギーの計算に用いられる. 以上のような各条件を LLG シミ
ュレータで設定して, 磁化解析をした.

3. 結果
解析は, サイズが 400×20×20[nm3]の直方体クラスタを想定して行った. 材料は, 鉄
鋼材料に用いられる鉄を想定して, 計算に使用する物性パラメータは
𝑀𝑠 =1714[emu/cm3] ， 𝐴 =1.5[μerg/cm][3], 𝐾 =4.7×105 [erg/cm3][4] と し た . void 欠 陥 は
2 × 2 × 2[nm3 ]のサイズで, クラスタ左端から100[nm]離れた所の中央部に, 以下のFig.1 に 示 す よ う に 配置 し た . 外 部 磁 場 は , 印 加 時 間 を 0 ～ 5600[ps], 大 き さ を-10000[Oe]～10000[Oe]の間で変化させて Fig.2 に示すように, クラスタの長手方向に
印加した.

Fig.1

Cluster with a void defect

Fig.2

Applied external field

また, 境界条件は, 磁壁の発生を 1 つにして解析しやすくするために, クラスタの両端の磁化を外側に固定するように設定した. そして, 初期磁化はランダムになるように設定し, 計算を行った.

93

JSST 2021
結果として, 磁壁が欠陥を通過する 5148[ps]のときに, void 欠陥の有無による変化
が, 交換エネルギー𝐸𝑒𝑥 と磁気双極子エネルギー𝐸𝑑𝑝 にみられた. 磁壁が通過する時間付近の, void 欠陥の有無による𝐸𝑒𝑥 の変化を Fig.3 に, 𝐸𝑑𝑝 の変化を Fig.4 にそれぞれ示
す.冒頭の説では, 磁壁が void 欠陥を通過する際には𝐸𝑑𝑝 は減少するとしていたが, 逆に増加していて, 通過直後の 5154[ps]のときには下がっている.一方, 𝐸𝑒𝑥 は磁壁がvoid 欠陥を通る際には下がり, 通過直後にはあがっている. これはミクロスケールのクラスタや欠陥を想定して計算を行ったため, 𝐸𝑑𝑝 より量子力学的なエネルギーである𝐸𝑒𝑥 の方が磁壁のトラップに影響を及ぼしたと考えられる. よって今後は欠陥の大きさを大きくしたり, 欠陥の数を多くしたりした場合, 𝐸𝑒𝑥 と𝐸𝑑𝑝 にどのような影響がみられるかを解析していくことを考えている. 実際, 磁性体内の欠陥の数や大きさは様々であり, より多くの void 欠陥を想定して解析を行うことで, 現実的な磁性体を想定して非破壊検査に生かすことが期待できる.

Fig.3 Changes of 𝐸𝑒𝑥 with or without a defectFig.4 Change of 𝐸𝑑𝑝 with or without a defect

参考文献[1] 古屋泰文, バルクハウゼンノイズ法による非破壊材料評価, 溶接学会誌, 64 巻,
pp.120-125, 1995.
[2] Louis Néel, Effet des cavités et des inclusions sur le champ coercitif, Cahiers de physique,N°25, pp.21-44, 1944.
[3] 近角聰信, 強磁性体の物理(下), 裳華房, p.175, 1995.
[4] 近角聰信,太田恵造,安達健五,津屋昇,石川義和, 磁性体ハンドブック, 朝倉書店,
p.779, 2000.

94

JSST 2021

水素脆化を受けたオーステナイト系ステンレス鋼の機械的特性評価
Mechanical properties of hydrogen embrittled
austenitic stainless steel
内田圭祐*, 山口克彦
Keisuke Uchida*, Katsuhiko Yamaguchi

福島大学 共生システム理工学研究科
Faculty of Symbiotic Systems Science,Fukushima University
s2170008@ipc.fukushima-u.ac.jp

Abstract.In this study, we performed a finite element analysis of hydrogen embrittlement. The value of the nucleation void volume ratio in the Gurson model was modified
to represent hydrogen embrittlement. To evaluate the analytical results, the lengthvariation of the width of the parallel section of the specimen and the fracture elongation were calculated. The length of the width of the parallel part of the specimen
showed a similar trend to that of the hydrogen embrittlement. The fracture elongation
also showed a similar trend to that of the hydrogen embrittlement case. The results
suggest that the hydrogen embrittlement test can be qualitatively described by varying the nucleation void volume ratio.

Keywords: Hydrogen embrittlement, Austenitic stainless steel,FEM,Gurson

1. 諸言
水素ステーションの実現に向けてさまざまな研究が行われている．その中で，材料に水素が侵入することで強度と延性を低下させる水素脆化という現象が問題であり，この現象に対する材料の強度の研究が行われている．水素脆化が顕著であるほど，引張試験の前95

JSST 2021
後での断面積の変化量が増加しにくくなると考えられる．また，水素脆化が生じると破断ひずみが低下するという研究報告がある 1)．水素脆化のメカニズムの一つとして水素助
長歪誘起空孔理論が提案されている．この理論は水素が侵入することにより転位の易動
度を増加させることでボイドの生成とその凝集を助長するというモデルである 2)．先行研
究では，この理論に基づき低炭素パーライト・フェライト鋼の 3 点曲げ試験による R 曲
線の水素脆化の有限要素法解析が行われた 3)．実際の解析では，水素脆化を表現するためにボイドの生成・成長を表す Gurosn モデル 4)を適用した 3)．本研究では，Gurson モデ
ルを用いて引張試験に関する水素脆化試験をオーステナイト系ステンレス鋼(以下，SUS304)で有限要素法解析を行えるか検討した．

2. 水素脆化試験解析
本研究で使用した有限要素解析のソフトウエアは MSC
Marc(MSC Software 製)である．このソフトではボイドの生成・成長を表す Gurson モデルが導入されている．水素脆化試
験の 2 次元解析モデルを図 1 に示す．図 1 は引張試験の JIS 規
格の 14B 号試験片の平行部である．境界条件は平行部の上端を拘束し下端を強制変位させるように設定した．本研究で使用した SUS304 の縦弾性係数 E は 195.85 GPa，ポアソン比 ν は0.29，降伏応力 σy は 332 MPa とする 5,6)．また， Gurson
モデルのパラメータは SUS304 の Gurson モデルに関する論文 図 1 解析モデル
を参照してパラメータを決定した．このときのパラメータを表 1 に示す 5)．q1，q2 は降伏曲面係数，f0 は初期ボイド体積比，fc は臨界ボイド体積比，fF は破壊ボイド体積比，εN は核生成平均ひずみ，SN は標準偏差を表している．本研究では，水素を侵入させたと仮定するためにボイド核生成の体積比 fN を 0.01～0.03 へそれぞれ変更した．
表1

Gurson モデルのパラメータ

q1

q2

f0

fc

fF

εN

SN

1.5

1

0.002

0.03

0.156

0.54

0.09

3.解析結果解析の評価を行うために，試験片の破断面の長さ ΔA と破断伸びλmax を算出した．ΔA
は次式(1)のように算出される．元の試験片の平行部の幅の長さ l，破断後の試験片の平行部の幅の長さを A とする．
∆𝐴 = 𝑙 − 𝐴

96

(1)

JSST 2021
式(1)から算出した結果を表 2 に示す．試験片の標点間距離を 25.6mm とし，水素脆化前
後の標点間距離の変位から算出した破断伸びλmax を表 3 に示す．L は破断の際の標点間距離の伸び，λmax は破断伸びを表している．
表 2 破断面の長さを比較

0.01
3.45
2.55

fN

A [mm]
ΔA [mm]

0.02
3.58
2.42

0.03
3.77
2.23

表 3 破断伸びの比較

fN

L [mm]
λmax [mm]

0.01
35.9
10.3

0.02
35.0
9.4

0.03
34.8
9.2

表 2 から，ボイド核生成の体積比 fN の値が上がると破断面の長さ ΔA の値が小さくなることがわかった．これはボイド核生成の体積比 fN が増加したことで，ボイド体積比がより早く破壊ボイド体積比 fF まで達したからだと考えられる．また，試験片の厚み方向の破断面の長さの変化量は試験が一方向にしか力がかからない試験であるうえ，ポアソ
ン比νが一定であるので，ボイド核生成の体積比 fN が大きくなるほど小さくなると考えられる．よって，ボイド核生成の体積比 fN が大きくなるほど，試験前後の断面積の変化量は小さくなることが見込まれる．このことから，Gurson モデルのボイド核生成の体積
比 fN を変更すれば水素脆化を表現することができると考えられる．
表 3 から，ボイド核生成の体積比 fN が増加すると破断伸びλmax が低下することがわかった．これはボイド核生成の体積比 fN が増加したことで，ボイド体積比がより早く破壊
ボイド体積比 fF まで達したためだと考えられる．この結果から，ボイド核生成の体積比

fN を使用すれば水素脆化で生じる現象の一つである破断伸びλmax の低下を表現することができる可能性があると考えられる．

4.結言
本研究では，水素脆化試験を想定した解析を行うために Gurson モデルを用いて有限要素法解析を行った．解析結果から，Guroson モデルで水素脆化の絞り，破断伸びを表現することができる可能性があることがわかった．今後は，引張強さを比較するために荷重-伸び線図を算出したいと考えている．

97

JSST 2021

参考文献
1)

小出賢一，南孝男，安樂敏朗，岩瀬明宏，井上博之，250℃の高圧水素ガス中でのSUS304 鋼の水素脆化感受性，材料と環境，2014，523-527．

2)

桃谷裕二，低炭素マルテンサイト鋼の水素脆性破壊挙動におよぼす変形条件の影
響，京都大学，博士論文，2017．

3)

南雲道彦，水素が関与する破壊の特徴，Vol. 56，No. 4，2007．

4)

Hiroshi YOSHIDA，Michihiko NAGUMO，FEM Analysis of Crack Growth in Fracture
Transition Region for Steels with Different Void Nucleation Frequency.ISIJ International，Vol. 38，No. 2，1998，196-202．

5)

Jianpeng Dong，Shilong Wang，Jie Zhou，Chi Ma，Sibao Wang，Bo Yang，Experimental
and numerial investigation on the shearing process of stainless steel thin-walled tubes inthe spent fuel reprocessing，2019．

6)新潟県工業技術総合研究所，http://www.iri.pref.niigata.jp/topics/H29/29kin24.html，
21th，June，2021.

98

JSST 2021

NaI(Tl)シンチレーション検出器を用いた放射線源
分布測定のシミュレーション
Simulation for estimation of radiation source distribution measured by
NaI(Tl) scintillation detector
菅 蒼太, 山崎 はるか, 山口 克彦
Sota Suga, Haruka Yamazaki, Katsuhiko Yamaguchi1

福島大学 共生システム理工学研究科

Graduate School of Symbiotic Systems Science and Technology , Fukushima University*

s2170025@ipc.fukushima-u.ac.jp

Abstract. The purpose of this study is to safely store substances with high concentrationsof radioactivity generated b y the accident at the Fukushima Daiichi Nuclear Disaster.
In this study, we simulated the radiation field for a measurement system using four NaI(Tl) scintillation detectors. As a result, it clarified that how the spectral pattern changes
with and without the shield.
Keywords: Safety store, Radiation source distribution

1. まえがき
福島第一原子力発電所事故により，放射性物質である燃料と原子炉の構造体などが溶融後，冷却・固化高濃度の放射性物質を含む燃料デブリ(以下，燃料デブリ)が大量に生成された．近年，この燃料デブリを原子炉建屋内から安全に取り出し，保管するための容器の具体的な提案も行われてきている．燃料デブリから放出される放射線量を定量的に評価するためには，容器内で不均一に存在している放射性燃料と遮蔽材料の分布を明らかにすることが重要となる．
そこで，NaI(Tl)シンチレーション検出器(以下，NaI スペクトロメータ)を用いて，外部から放射線を測定する装置を試作している．この装置を使用して様々な方向からスペクトルを測定し，その強度に差が生じれば遮蔽材料や高線量を示す位置が推定でき，収納缶の内部構造も推定ができると予想される．今回のシミュレーションでは，複数本の NaI スペクトロメータを様々な方向に配置し，測定することで遮蔽材料の効
果の有無を判別すること、またスペクトルの測定に最適な本数や配置を決定することを目的としている。
以上を踏まえ，放射性燃料と遮蔽材料が不均一に分布したモデル試料を想定し，試料から放出されるγ線スペクトルの遮蔽材料における依存性をモンテカルロ計算コ99

JSST 2021
ード PHITS (particle and Heavy Ion transport code System)を用いて測定に先駆けて実施したシミュレーション結果を報告する[1]．

2. 方法
解析にあたり，放射線源は 137Cs(全吸収ピークのエネルギー：0.662MeV)を，遮蔽材
料には鉛(特性 X 線のエネルギー：0.083MeV)を用いた．
測定系を真上から見たものの一例を Fig.1 に示す。Fig.1 に示した NaI スペクトロメータは 4 本であるが、実際には 4 本以外に複数本使用している．二次元平面上の中心に 137Cs 線源を置いて，複数本の NaI スペクトロメータの検出部の先端から線源までの距離が一定になるよう配置した．また中心に置いた 137Cs の一部を遮蔽材料で覆い，複数のうちいくつかの NaI スペクトロメータを遮るように遮蔽材料を配置した。また
NaI スペクトロメータの本数や，線源に対する角度を変化させ，各方向からγ線スペ
クトルを測定した．二次元平面上でエネルギースペクトルを測定することにより，放射線源と遮蔽材料との関連性を検討した．

Fig.1 Measurement system using NaI (Tl) scintillation detector

3. 結果
Fig.2(a)に放射線源に対して遮蔽材料を用いていないもの，Fig.2(b)に放射線源に対して遮蔽材料を用いたものをそれぞれを示す．今回のシミュレーションで，137Cs の0.66MeV の全吸収ピークでは約 2 倍の変化が見られ，0.1MeV の低エネルギー側では，遮蔽材料が無い場合に現れていた 137Cs の 0.032MeV の特性 X 線ピークが，遮蔽材料がある場合には現れず、代わりに鉛の 0.083MeV の特性 X 線ピークが出現している。
下図からは鉛の遮蔽材料による変化が見受けられ，鉛の遮蔽材料による遮蔽効果は大きな影響を及ぼすことが分かった．今後，複数の NaI スペクトロメータから検出されるエネルギー強度の相対関係から，スペクトルパターンを解析していくことで高濃度のデブリを含む格納容器の内部状況が分かると予想される．

100

JSST 2021

Fig.2(a) gamma ray spectrum

Fig.2(b) gamma ray spectrum

not using a shieldusing a shield

4. 参考文献[1] Sato T, Iwamoto Y, Hashimoto S, Ogawa T, Furuta T, Abe S, Kai T, TsaiP, Matsuda N,
Iwase H, Shigyo N, Sihver L and Nitta K, “Features of Particle and Heavy Ion Transport
code System (PHITS) version 3.02,” Journal of Nuclear Science and Technology, 2018,56.

101

JSST 2021

構造材料を対象とした熱ルミネッセンス
測定のためのシミュレーション
Simulation of thermoluminescencemeasurement of structural materials
上杉 美咲*, 高瀬 つぎ子, 山口 克彦
Misaki Uesugi*, Tsugiko Takase, Katsuhiko Yamagushi

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

Abstract. A Monte Carlo simulation was performed to analyze the glow curve for calculating
the integrated radiation dose of structural materials using thermoluminescence. In our experiments, we found that the shape of the glow curve changed when the sample was irradiated fora long time. This study enabled us to calculate the glow curve of a thermo-fluorescent material
with a single trap. The glow curve of a thermo-fluorescent material with a single trap can becalculated by using this method.
Keywords: Boltzman distribution, Dosimeter, Thermoluminescence, Monte Carlo simulation

1. まえがき
外部から放射線などの刺激を受けた物質を加熱すると発光する熱ルミネッセンス(Thermoluminescence:TL)という現象がある．TL の発光量は放射線吸収量と比例する
ため，積算放射線量を推定することが可能であり，広島・長崎に投下された原爆から出たγ線量[1]の算出にも用いられてきた．2011 年 3 月に発生した福島原子力発電所の事故においても，事故発生から 10 年が経過しており，当時の空間線量を直接求めることが困難になっている．そこで，発電所周辺の構造材料から熱蛍光物質を採取し TL法を用いることで，原爆同様，放射線量を間接的に推定可能であると考えた．
積算放射線量は，グロー曲線と呼ばれる熱蛍光物質を加熱したときの熱蛍光物質温度とその温度における熱蛍光量の関係[2]より算出されるため，測定系の精度が要求される．製作した測定系の精度を調べるため，一般的な熱蛍光物質である CaSO4:Tm の
TL 測定を行った．CaSO4:Tm は，β線とγ線の年間線量決定に高感度で安定性が良い素子であるという特徴[3]をもつため，TL の測定に多く用いられている．

102

JSST 2021
CaSO4:Tm の TL 測定より，放射線照射直後に測定した試料と，照射後 5 ヶ月経過したのちに測定した試料を比較すると，グロー曲線に違いが見られた．これを Fig.1に示す．照射直後に測定した試料においては２つのピークが確認できたものの，照射後 5 ヶ月経過したのちに測定した試料では 1 つのピークしか確認できなかった．時間経過により熱蛍光物質のトラップ準位が変化したと考えられる．照射後の時間経過によるグロー曲線への影響が分かれば，照射後の保管期間に応じて適切な解析を行うことができる．本研究では，積算放射線量推定の精度確立のために，試料照射後の時間経過によるグロー曲線の変化を解析することを目的とした．

Fig.1 Glow curve of CaSO4:Tm

2. 手法
Fig.2 にエネルギーバンドモデルによる熱蛍光過程を示す．熱蛍光物質に放射線が当たると，相互作用により電離する．電離した電子は励起され価電子帯から伝導帯へと移動する．伝導帯の電子は捕獲中心に捕捉され準安定状態となり，ホールは物質中に捕捉され蛍光中心となる．
この状態の物質を加熱すると，電子が捕獲中心から伝導帯へと移動し，蛍光中心のホールと再結合する．このとき，伝導帯と蛍光中心のエネル
ギー差が TL として発光する．
本研究では，捕獲中心から伝導帯へと移動した電子の数を TL 発光強度とし，ボルツマン分布を用いて算出した．ボルツマン分布を(1)式に示す．

𝛪(𝑇) = 𝑛! exp (− 𝐸 ⁄𝑘𝑇)

Conduction band

electron trap
heat

TL

hole trap

Valence bandradiation

Fig.2 Energy band model
scheme of TL(1)

ここで，I (T)は温度 T での TL 強度，n0 はトラップされている総電子数，E は伝導帯と蛍光中心のエネルギー差，k はボルツマン定数である．この式より確率条件を設定し，モンテカルロ法を用いて TL の発光機構を解析した．

103

JSST 2021

3. 結果と考察
伝導帯と蛍光中心のエネルギー差 E が 0.10~0.15eV の場合のグロー曲線推定結果を Fig.3 に示す．TL の発光機構の通り，エネルギー準位差に依存してピーク温度が異なることが確認でき，1 つのトラップを持つ熱蛍光物質のグロー曲線を算出することができた．この推定結果を用いると，CaSO4:Tm は 0.11eV と 0.14eV のトラップをもっており，照射後の時間経過により 0.11eV のトラップ準位が移動したと考えられる．この手法を応用すると，CaSO4:Tm のように２つのトラップを持つ熱蛍光物質のグロー曲線の算出ができ，時間経過によるトラップ準位の移動過程の解析を行えると期待できる．

Fig.3 Glow curve by simulation

参考文献[1] 市川米太,熱ルミネッセンスによる原爆線量の測定 I 広島•長崎の被爆瓦の示す発光曲線,奈良学芸大学紀要.自然科化学,12,pp.51-57(1964).
[2] 眞 正 浄 光 , 古 場 裕 介 , 熱 蛍 光 線 量 計 の 諸 特 性 と 応 用 研 究 の 紹 介 , 放 射 線 科
学,103,p.13(2017).
[3] 長友恒人,山本健司，熱ルミネッセンス年代測定法におけるガンマ線年間線量測定 法 の 簡 便 化 - 市 販 TLD カ プ セ ル を 用 い る 試 み -, 古 文 化 財 教 育 研 究 報
告,15,pp.31-39(1986).

104

JSST 2021

非定常高周波電磁界解析手法の開発と応用の検討

Study for Development and application of Nonsteady
High-Frequency Electromagnetic Field method

西村 聡汰 1*
Souta Nishimura1*
武居 周 1
Amane Takei1

11

宮崎大学 工学部 電気システム工学科

Department of Electorical and Systems Engineering, Faculty of Engineering,Miyazaki University*

hl18034@student.miyazaki-u.ac.jp

Abstract. The high-frequency electromagnetic fields solver: ADVENTURE_FullWave is a
finite element analysis software in the parallel computational mechanics system: ADVENTURE System. The ADVENTURE_FullWave can perform detailed, high-speed, andsteady-state high-efficiency finite element analyses of large-scale electromagnetic field problems using the hierarchical domain decomposition method (HDDM) and the corresponding
parallel distributed processing environment. However, also non-steady analysis is required. In
this research, we study and develop the non-steady high-frequency electromagnetic field
method. The detail of non-steady high-frequency electromagnetic field analysis will be discussed on the patch antenna problem in the JSST2021 conference.
Keywords: Finite element method, Electromagnetic fields, Large-scale analysis, Hierarchicaldomain decomposition method

1. 背景
近年の計算機の性能向上と数値計算技術の進歩により，計算機を用いた，高周波帯域の電磁界解析が普及した．これにより，アンテナ等の開発・ 設計において，実際に機器を試作し，測定を行わなくとも，電磁界の振る舞いを予測することが出来るよ
うになり，開発期間の短縮や開発コストの削減に非常に有用なツールとして広く利用されるようになった．オープンソースソフトウエア ADVENTURE_FullWave[1]もその１つであり，解析モデルの作成・メッシュ生成・境界条件設定などで ADVENTURE

105

JSST 2021
System の各種モジュール群を基盤に用いた解析システムの構築とともに運用される．
本解析コードは，高精細モデルを用いた定常問題の大規模解析が行える利点がある
一方，非定常計算に対応していないため，電子回路設計等で EMC などの評価に適用することができない．
そこで，本研究ではまず有限要素法に基づく高周波電磁界解析の非定常計算機能について検討し，今後の開発指針を得る．

2. 基礎式
2.1. 時間変動場のベクトル波動方程式
電磁界の基礎方程式である Maxwell の方程式のうち，(1a)式の Faraday の法則，
(1b)の Ampere−Maxwell の法則と，電界，磁界に関連する構成則(1c)により時間変動
場Ωの電界𝑬を未知関数としたベクトル波動方程式(2)が導かれる[1]．
𝑟𝑜𝑡 𝑬 = −

𝜕𝑩
𝜕𝑡

(1a)

𝜕𝑫
𝜕𝑡

(1b)

𝑫 = 𝜀𝑬, 𝑩 = 𝜇𝑯

(1c)𝑟𝑜𝑡 𝑯 = 𝑱 +𝜕𝑬𝜕𝑱
𝑟𝑜𝑡 2315𝜇6𝑟𝑜𝑡 𝑬7 + 𝜀=−𝜕𝑡𝜕𝑡
𝑬 ×𝒏=0

in Ω

(2a)
(2b)

on ∂Ω

@
𝑱 = σ𝑬(2c)

ここで，𝑬は電界, 𝑯は磁界, 𝑫は電束密度, 𝑩は磁束密度, 𝜀は誘電率, 𝜇は透磁率, 𝑱は電流密度, 𝒏は境界における外向き法線ベクトル, 𝑡は時刻である．解くべき方程式は式(2a)であり，電界𝑬を未知数とした E 法として定式化される．本定式化において，透磁率𝜇，複素誘電率𝜀を既知数とし，放射源として電流密度𝑱を直接，あるいは式(1c)@ を与える．また，得られた電界𝑬より(1a)式を用いて磁界𝑯を計算する．
より既知電界𝑬2.2. 有限要素定式化
式(1a)より弱形式を導き，領域Ωの有限要素分割を考える．電界 E を Nedelec 四面体一次要素（辺要素）で近似し，電流密度 J を通常の四面体一次要素で近似する．
𝐄! , 𝐉! をそれぞれ E , J の有限要素近似とする．𝑬∗" は任意の試験関数とする[2]．
! (1⁄𝜇 ) rot 𝐄! ∙ rot 𝐄!∗ dv +Ω𝜕2
𝜕! ε𝐄! ∙ 𝐄!∗ dv = −! 𝐉! ∙ 𝐄!∗ dv
2
𝜕𝑡
𝜕𝑡
Ω
Ω

(3)

式(3)を解くべき有限要素方程式とする．時間発展は Newmark-b法による陰的解法に

106

JSST 2021より扱う[3]．

3. 数値例
性能評価に用いる計算機は，筆者らの研究室で所有するマルチコア PC である．本
計算機は，Intel Core i7-2600K のマルチコア CPU および 32GB のメモリが搭載されている PC であり，並列化ライブラリとして MPI(Message Passing Interface)を用いる．図 1(a)に示す要素数 10 万程度の円形パッチアンテナを波源として，(3)式に基づく電磁界解析を行った結果を図 1(b)に示す. 本数値例では，パッチアンテナに対して図1(b)に示す電圧波形を得た．本数値例に関して式(3)を導入して解析を行った。詳細な議論は講演にて行う．

(a) 円形パッチアンテナの解析(b) 得られた電圧波形

図 1 円形パッチアンテナの解析

参考文献[1] Jian-Ming Jin: The Finite Element Method in Electromagnetics – Third edition, Willey,
2014[2] A. Takei, S. Sugimoto, M. Ogino, S. Yoshimura and H. Kanayama: Full Wave Analysesof Electromagnetic Fields with an Iterative Domain Decomposition Method, IEEE
Transactions on Magnetics, Vol. 46, No. 8, pp.2860-2863(2010) DOI:
10.1109/TMAG.2010.2044775.
[3] 村上裕哉，山本広太，工藤彰洋，武居周, “時間領域並列有限要素法に基づく大規模音響解析”,日本シミュレーション学会論文誌, Vol.10, No.4, pp.89- 98 (2018)
DOI: 10.11308/tjsst.10.89.

107

JSST 2021

混雑度と境界条件による大小粒子の局在化現象：
cm スケールでの細胞モデル実験
Localization phenomenon of large and small particles
by congestion degree and boundary conditions :Cell model experiment on cm scale
黒田真帆 1*, 鷹取慧 1, Chwen-Yang Shew2, 剣持貴弘 1, 吉川研一 1
Maho Kuroda1*, Satoshi Takatori1, Chwen-Yang Shew2, Takahiro Kenmotsu1 Kenichi
Yoshikawa3

1

同志社大学 生命医科学研究科
2
ニューヨーク市立大学
1
Facul. of Life and Med. Sci. Doshisha Univ.
2
Dept. of Chem. City Univ. of New York
*ctuf0018@mail4.doshisha.ac.jp

Abstract. The various molecules and organelles in living cells are suitably positioned within
the cell to carry out their functions at the appropriate time. The nucleus moves toward andmaintains its position at the center of the cell. We conducted a real-space experiment inwhich the intracellular structure was modeled by large and small particles. In our experiments,
we focus on that the inside of the cell is very crowded and that there are non-equilibrium
fluctuations. From the results of this experiment, it was clarified that the behavior of particles
differs depending on the degree of congestion and the softness of the boundary conditions.
Keywords: Self-organization, Non-equilibrium fluctuation, Molecular dynamics simulation

1. 序論
細胞内の生体分子は秩序構造を維持することで正常に機能している．細胞内の秩序形成メカニズムを明らかにすることは，細胞を自律的に機能する生命システムとして捉えることを可能とし，生命科学だけでなく，医学・医療分野の基礎的な知見となることが期待される．しかし，細胞内の生体分子の自律的な秩序形成メカニズムについては未解明な点が多い．

108

JSST 2021
細胞内秩序構造の中でも本研究で注目しているのが，細胞核の細胞内配置である．
多くの真核細胞において細胞核は細胞の中央付近に位置している[1]．また，細胞内は，生体高分子が 1g 中に 0.3～0.5g 存在する混雑した環境であるということと，常にエネルギーの消費が行われており，著しい非平衡ゆらぎに晒されていることが細胞内の生体高分子の秩序だった構造や機能に直接関係していると考えられている[2]．本
研究では，実際に細胞内の生体高分子の挙動を観察することは困難であるため，細胞内の環境を模擬できる，加振器を基軸とした実空間モデルを採用する．本研究では，生体高分子に見立てた大きさの異なる多数の粒子が，ゆらぎを受けながら，どのように自発的に秩序を形成するかを明らかにする．また，分子動力学やモンテカルロシミュレーションを実施し，理論的な考察も行い，実験結果との検証から物理モデルを構築する．

2. 実験方法
実験には主に大球，小球，お皿のような容器，加振器を用いる．壁面と底面が曲面で繋がれた容器を使用し，容器内に直径 15mm の大球 1 つと直径 6mm の複数の小球を入れ，加振器を用いて振動を加える．加振は周波数を 60Hz，振幅を 0.25mm，
最大加速度を 18m/s2 とする．
図 1 に加振器の実験系と実験条件をまとめた図を示す．

図 1 (a)加振器の実験系と加振のデータ (b)容器，小球の寸法
このように容器の壁面と底面が曲面で繋がれていることにより，柔軟性によるバネのような効果が現れることを期待した．
容器底面部に小球を敷き詰めた時（小球が 150 個）を混雑度 100％に設定し，小球の数を変えることで混雑度 20%，40%，60%，80%の状態を作る．また，加振はいず
れも 300 秒間行い，0 秒の時は大球を小球のクラスターの端部に置く．

3. 研究結果
加振時に大球が示す軌跡を時間的に可視化したものを図 2 に示す．図 2 から混雑度が高くなるにつれて大球の軌跡がフラットな底面と曲面の境界である黒点線付近に集中していることから，混雑度が高いほど大球は境界付近に局在化していることがわ
かる．
また，内部と外部での存在比を order parameter とし解析したところ，今回の実験系
では，臨界点近傍の振る舞いが現れることも明らかとなった．その図を図 3 に示す．
図の横軸は容器中心からの距離を表しており，左端が容器中心である．図 3 から，混
雑率 20%，40%では中心側でグラフが落ち込んでいるが，混雑率 60%で臨界近傍の様子を示し，混雑率 80%で外側でグラフが落ち込んでいる．

109

JSST 2021

図 2 各混雑度における大球の軌跡の時間推移。
黒実線が容器の壁面、黒点線が容器底面の境界を表す。

図 3 大球の自由エネルギー曲線

4. まとめ本研究では，大小粒子集団にゆらぎを印加した場合に混雑度と境界条件によって大球の局在位置がスイッチングすることを見出だしている．その結果は Shew らによるモンテカルロ法のシミュレーションにより得られた結果と同様になったと考えることができる[3]．また，本研究では cm スケールでの実験を行っているので、加振器を用いてゆらぎを発生させており、細胞スケールでの Brown 運動とは、ゆらぎの起源が異なっている。このような違いにも拘わらず、スケールの異なる実験系が同様の傾向を示していることは興味深い。本研究を実施することで得られる知見は，細胞を構成する生体高分子がゆらぎを受けながらどのように自発的に秩序を形成し，生体機能

110JSST 2021を発揮しているのかについてその仕組みを明らかにし，生命現象の本質的な理解を促進するものである．

参考文献[1] Tanimoto, H., Kimura, A., Minc, N., Shape–motion relationships of centering
microtubule asters, Journal of Cell Biology, vol.212, no.7, (2016), pp.777-787.
[2] Oda, S., Kubo, Y., Shew, C. Y., & Yoshikawa, K. (2016). Fluctuations induced transitionof localization of granular objects caused by degrees of crowding. Physica D: Nonlinear
Phenomena, 336, 39-46.
[3] Shew, C. Y., Oda, S., & Yoshikawa, K. (2017). Localization switching of a large object
in a crowded cavity: A rigid/soft object prefers surface/inner positioning. The Journal ofchemical physics, 147(20), 204901.

111

JSST 2021

Sightseeing Guidance System to Maximize
Satisfaction Using Real-Time Spot Information∗
Yuya Sato1 , Yoko Nakajima2 , Hirotoshi Honma2,*
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∗ honma@kushiro-ct.ac.jp

Abstract. In this study, we challenged the provision of sightseeing routes using real-time
information such as weather and congestion. This is a system that allows tourists to input
their wishes and provides them with the best sightseeing route based on their wishes andreal-time information, using graphs and optimization problems. We succeeded in providing
a tourist route using this method, but confirmed that it is still far from practical use.
Keywords: Graph theory, Discrete optimization, Operations research

1.

Introduction

The Eastern Hokkaido region, in which Kushiro City is located, is an area of impressive
natural beauty, which features the Kushiro Marshes, the largest marshlands in Japan, as
well as both the Akan and Mashu national parks. Additionally, it is home to the historicallyrecognized Ainu culture, which has received a significant amount of attention in recent years
from many quarters. The region has received recognition for the national government’s
measures in regard to the Mizu-no-Kamui tourist area and Wide-area tourist route, which
promote tourism.
Currently, Kushiro City provides tourists with the Kushiro Akan Mashu Hospitality Navi
APP., which enables tourists to use LPQA and free public Wi-fi, regional bus company
shared locations, and AI chat bot services focusing on tourism. However, a current issue isthat data with diﬀerent properties are being used in independent environments, and there isno collaboration with, or utilization of, other services.
In this study, we imagine a virtual environment in the near future, in which a high-speed
communication network is provided in the Eastern Hokkaido region, thus enabling weather,temperature, humidity, specific event, and congestion information for each sightseeing spot∗

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.

112to be acquired in real time, and we undertake the challenge of developing a sightseeing
guidance system that can maximize the individual satisfaction levels of each tourist visiting
Eastern Hokkaido under this environment.

2.

System Overview

This system consists of local information on the weather, temperature, humidity, specific
events (content/time/cost), and congestion, etc., for the various sightseeing spots, a databasefor setting the satisfaction levels of travelers, which fluctuate in combination with thesefactors, and map data that include the route and necessary time to move between each sightseeing spot.
By inputting the requirements of the user into the system (departure/arrival location, departure/arrival time, time allowed for sightseeing, meal requests, costs), it provides a sightseeing plan that will maximize the satisfaction level of the user within the time allotted for
tourism. An image of the system is shown in Figure 1.

Figure 1: Image of route provision
A mathematical optimization model was used to derive the optimal tourist route. This
reduces the confusion experienced by the traveler in selecting from multiple routes provided
by guidebooks. Additionally, as, depending on the sightseeing spot, there may be cases inwhich the satisfaction level fluctuates based on the weather, level of congestion, and events,
in this way, the optimal route can be provided, and the traveler can enjoy sightseeing.

3.

Current Situation and Future Outlook

To date, we have completed the creation of an optimization model for providing sightseeing
routes and a database of sightseeing spots in Kushiro City and other areas in the eastern
Hokkaido region, which consists of basic information about these spots. Based on these
data, we have also succeeded in providing tourist routes for each request.
Future tasks include the development of an interface and system for actual use. In addition, it is necessary to acquire and operate real time data to the extent possible. The final
goal is to make the information of tourist attractions more realistic by having actual tourists
use the system.

113

JSST 2021

DNA のトリチウム耐性検討のための DNA テロメ
ア構造中水素の溶媒接触表面積の計算
Evaluation of solvent-accessible surface area of backbone hydrogen in telomeric DNA for studying tritium
resistance of DNA
土田 陽平 1*, 齋藤 誠紀 1, 中村 浩章 2, 米谷 佳晃 3, 藤原 進 4
Yohei Tsuchida1*, Seiki Saito1, Hiroaki Nakamura2,
Yoshiteru Yonetani3, Susumu Fujiwara41

山形大学
核融合科学研究所
3
量子科学技術研究開発機構
4
京都工芸繊維大学
2

Yamagata University
National Institute for Fusion Science3
National Institutes for Quantum and Radiological Science and Technology
4
Kyoto Institute of Technology
12*

t212989m@st.yamagata-u.ac.jp

Abstract. Tritiated water is generated under the decommissioning process of the Fukushima
Daiichi Nuclear Power Station. In addition, tritium is planning to use as fuel in fusion power
plants, which is expected as a future power generation technology. Therefore, it is important tounderstand the impact of tritium on biomolecules in living organisms including human in detail.
We aim to elucidate the mechanism of DNA damage due to the radioactive decay effect that
occurs when light hydrogen in human DNA is replaced with tritium, using molecular dynamics(MD) methods. To understand the decay effect on DNA, first, it is necessary to evaluate thereplaceability of light hydrogen to tritium for each hydrogen in DNA. In this study, to evaluate the degree of replaceability of the backbone hydrogen atoms in telomeric DNA of human,
solvent-accessible surface area (SASA) is calculated for the data obtained by MD simulation.
As a result, it is found that the SASA of H5 hydrogen is large in the hydrogen atoms in the
backbone.
Keywords: DNA, solvent-accessible surface area, tritium, radioactive decay effect

1. 研究背景・目的
現在、東京電力福島第一原子力発電所で発生したトリチウムの残留している処理水

114

JSST 2021の処分方法について、多くの議論が交わされている。また、核融合発電は、重水素とトリチウムの核融合反応を用いた発電方法である。燃料として使用される重水素とリチウムは海水から無尽蔵に採取できる。また発電の際、高レベルの放射性廃棄物が発生しないため、放射線被ばくの危険性も低い。このような事情から、トリチウムの生体分子への影響を詳細に解明することが求められている。本研究では、ヒト DNA 中の軽水素がトリチウムに置換した際に生じる壊変作用による DNA 損傷のメカニズムを、分子動力学法を用いて解明することを目指している。そこで、分子動力学法を実施する前段階として、トリチウム置換しやすい場所を特定するため溶媒接触表面積を計算した。

2. シミュレーションモデル [1]
水分子で満たした100Å×100Å×100Åの計算領域に、ヒトDNAのテロメア構造を配置する。テロメア構造の配列はTCTAGGGTTAGGGTTAGである。DNAの電荷を中和するため、イオン濃度が0.15 mol（生体内の濃度）となるようにNa+とCl−を追加しイオン化を行った。シミュレーションモデルを図１に示す。

図1

DNAシミュレーションモデル[1]周囲の水分子は非表示とした。MDシミュレーションには、LAMMPSを用い、力場にはReax-FF(CHON-2017_weak)を用いる。x, y, z方向には周期境界条件を課した。粒子数
N、圧力P、温度Tが一定の条件下（NPTアンサンブル）で計算した。熱・圧力浴はNoséHooverを使用した。人体中の温度・圧力に近くなるように、設定圧力、温度はそれぞれ1.0 atm、310 Kとした。時間ステップは0.25 fsとし、100 psの間の運動を計算した。
ここで、水分子で囲まれた DNA の水素における OH ラジカルの反応率と溶媒中の水素における溶媒接触表面積には相関関係があることが知られている[2]。本研究でも、ヒト DNA 中の溶媒接触表面積がトリチウム置換のしやすさと相関があると考え、
ヒト DNA 中バックボーン水素について MD 計算から得た原子の位置座標を用いて溶媒接触表面積を計算した。

3. 計算結果 [1]溶媒接触表面積の時間平均を、DNAの2本鎖strand1、strand2ごとに計算したうちの

115

JSST 2021
strand2の計測結果を図2に示す。

図2 溶媒接触表面積の平均値[1]命名規則は文献[3]に従っており、H2およびH5の溶媒接触表面積は、H2'とH2''および H5'とH5''の溶媒接触表面積の値を足し合わせて求めた。エラーバーは標準偏差を示している。また、DNAバックボーン部におけるH5水素の溶媒接触表面積の時間平均を計算し、どの塩基においても、H5'の溶媒接触表面積が大きいとH5''の溶媒接触表
面積も大きくなる傾向があることが確認できた。
本予稿では、DNA構造の熱振動による溶媒接触表面積の動的変化について論じた。
ポスターセッションでは、DNA中水素の溶媒接触表面積とDNA中水素の配置の関係を解明することを目的に、任意の塩基配列のDNAの軽水素における溶媒接触表面積の動的変化について論じたい。

参考文献[1] 土田陽平, 齋藤誠紀, 中村浩章, 米谷佳晃, 藤原進: 分子動力学法を用いた
DNA テロメア構造中バックボーン水素の溶媒接触表面積の評価, 日本シミュレーション学会論文誌, 13, (2021), 32-26.
[2] Y. Yonetani, H. Nakagawa: Understanding water-mediated DNA damage production bymolecular dynamics calculation of solvent accessibility, Chem. Phys. Lett., 749 (2020),
137441.
[3] D. Sy, C. Savoye, M. Begusova, V. Michalik, M. Charlier, and M. Spotheim-Maurizot:
Sequence-dependent variations of DNA structure modulate radiation-induced strand
breakage, Int. J. Radiat. Biol., 72, (1997), 147-155.

116

JSST 2021

新規な病態計測法：組織の伸展応答パターンの活用

Smart Diagnosis on Disease State based on
Characteristic Cracking Pattern by Stretching
Response of Tissue Slice田上 幸歩 1,* , 鷹取 慧 1 , 剣持 貴弘 1 , 吉川 研一 1
Tagami Yukiho1,* , Satoshi Takatori1 , Takahiro Kenmotsu1 , Kenichi Yoshikawa1
1 同志社大学大学院 生命医科学研究科 医工学・医情報学専攻 医工学コース
1 Faculty of Life and Medical Sciences, Doshisha University∗ ctuf0033@maail4.doshisha.ac.jp

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

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

117