CODE2015
3rd INTERNATIONAL CONFERENCE
ON COMPUTATIONAL DESIGNIN ENGINEERING
December 14-15, 2015
Tokyo, Japan
Organized byComputational Structural Engineering Institute of Korea (COSEIK)&International Center for Urban Engineering (ICUS), IIS, UTokyo

Supported by
Institute of Industrial Science (IIS), the University of Tokyo (UTokyo)

Table of Contents

Welcome Address ···························································· 1
Organization of CODE2015 ················································ 2General Information on CODE2015 ······································· 3
Program Overview ··························································· 7
Keynotes······································································ 21
Technical Program – Day 1 Session 1·································· 27
Technical Program – Day 1 Session 2·································· 35
Technical Program – Day 1 Session 3·································· 43
Technical Program – Day 1 Session 4·································· 51
Technical Program – Day 2 Session 5·································· 57
Technical Program – Day 2 Session 6·································· 63
Technical Program – Day 2 Session 7·································· 69
Technical Program – Day 2 Session 8·································· 75
Technical Program – Day 2 Session 9·································· 81
Technical Program – Day 2 Session 10 ································ 87

Welcome Address
It is my great honour and pleasure to invite you to the 3rd International Conference on Computational
Design in Engineering (CODE2015) to be held from December 14th to 15th, 2015 at the Institute of Industrial Science,
The University of Tokyo, Japan, hosted by the International Center for Urban Safety Engineering (ICUS) andorganized by the Local Organizing Committee of CODE2015.
The remarkable development of Computer Science in recent years has contributed greatly to theadvancement of research in various fields including the field of engineering. Conducting an actual experiment or
survey could be difficult due to limitation in scale or time constraint, material constant, or boundary condition
however by utilizing computer simulation, understanding of phenomena and proposal of problem-solving measures
are making rapid progress.
Under the theme “Current and Future Directions in Engineering with Computers”, CODE2015 will focus
on the Computational Engineering and Mechanics and IT disaster management. CODE2015 will be an excellent
learning experience in the field of ‘Computational Design and Engineering’ for all participants. Globally renowned
speakers as well as many professionals and experts from around the world have been invited to CODE2015.
CODE2015 will be the third conference following the great success of CODE2009 and CODE2012. Join
us to learn and share the most recent advances in the related fields and improve your social network in our eversmaller global community. The Organizing Committee welcomes all the participants with Japanese Omotenashi(hospitality) spirit and hope to provide full support during the conference.
I am committed to making sure that this conference will be a valuable and enjoyable event. I would liketo express my deepest gratitude to all of you for attending CODE2015, to all sponsors for their generous support
ensuring the success of this conference, and lastly, to those who made this conference possible with their unwavering
dedication.
Once again, all of the members of CODE2015 Organizing Committee welcome you to CODE2015 andTokyo, Japan.

Prof. K. Meguro
Chairman
Organizing Committee of CODE2015―1―Organization of CODE2015
Co-Chairman
Waon-Ho Yi, Professor, Kwangwoon Univerity, Korea
Kimiro Meguro, Director, ICUS, IIS, UTokyo, Japan

Steering Committee
Heoung-Jae Chun, Professor, Yonsei University, Korea
Sang-Ho Lee, Professor, Yonsei University, Korea
Chee-Kyeong Kim, Professor, Dankook University, Korea
Maeng-Hyo Cho, Professor, Seoul National University, Korea
Hong Huang, Professor, Tsinghua University, China
Teruo Fujii, Director General, IIS, UTokyo, Japan
Yoshiaki Nakano, Professor, IIS, UTokyo, Japan
Chisachi Kato, Professor, IIS, UTokyo, Japan

Technical Committee
Hyung-Joon Kim, Associate Professor, University of Seoul, Korea
Hyuk-Chun Noh, Associate Professor, Sejong University, Korea
Dae-Jin Kim, Assistant Professor, Kyunghee University, Korea
Sang-Joon Shin, Professor, Seoul National University, Korea
Marie Oshima, Professor, IIS, Utokyo, Japan
Kazuo Kashiyama, Professor, Chuo University, Japan
Kohei Nagai, Associate Professor, ICUS, IIS, UTokyo, Japan

Organizing Committee
Muneyoshi Numada, Lecturer, ICUS, IIS, UTokyo, Japan
Ho Choi, Assistant Professor, IIS, UTokyo, Japan
Hideomi Gokon, Assistant Professor, IUCS, IIS, UTokyo, Japan

Contact
Eiko Yoshimoto―2―

GENERAL INFORMATIONLANGUAGE
The conference will be conducted in English.

SYMPOSIUM VENUE
Convention Hall, Komaba Research Campus
As311,312
Tel:

As313,314+81-3-5452-6472

As303,304
Fax:+81-3-5452-6476

Email: code2015@iis.u-tokyo.ac.jp
Web: http //icus.iis.u-tokyo.ac.jp/CODE2015.html

ACCESS to Komaba Research Campus―3―Map of IIS Campus―4―Inside map of conference venue―5―

REGISTRATION DESK
The registration desk will be open in front of the Convention Hall for all participants andaccompanying guests during the following times:
December 14

12:00 - 13:30

December 15

08:30 - 08:50

The registration fee for participants and students includes one volume of abstract and symposium
program, USB proceedings,1 lunch, 3 tea-breaks, and 1 banquet
Foreign Delegates:
USD300 (JPY30,000) Early Bird Registration (Before October 15, 2015)
USD350 (JPY35,000)
Students: (USD100) JPY10,000

WELCOME RECEPTION
The Conference Welcome Reception will be held at Capo Pericano at Komaba campus at 18:00 on
Tuesday, December 14, 2015. We will provide the coupon for the welcome reception in your conference
bag. Additional tickets for accompanying persons can be purchased from the Registration Desk for
USD 40/person.(JPY4500)

LUNCH
The Symposium Lunch will be held at Foyer of Convention Hall on December 14 and 15, 2015. We
will provide the coupons for lunch in your conference bag. Additional coupon for accompanying
persons can be purchased from the Registration Desk for USD20/person.(JPY2500)

INSTRUCTIONS FOR SPEAKERS
Keynote Speech: 20 minutes presentation
Common Speaker (Presentation):
15 minutes including 12 minutes presentation and 3 minutes discussion
A common computer will be available for speaker presentations if necessary. Because of computer
security reasons, we would appreciate it if speakers prepare their own computer.

―6―

3rd International Conference on Computational Design in Engineering

PROGRAM OVERVIEW
Time

Monday, 14 December

12:00-13:30

Registration (Foyer 2nd floor, An block)
Opening address (Convention Hall)
1. Welcome Speech (Prof. Teruo Fujii, Director of IIS, The University of Tokyo, Japan)

13:30-13:50

2. Welcome Speech (Prof. Waonho Yi, Dept. of Architectural Engineering, Kwangwoon
Univ., Seoul, Korea)
3. Opening Remark (Prof. Kimiro Meguro, Director of ICUS, IIS, The University ofTokyo, Japan)
Moderators: Dr. Hideomi Gokon
Keynote Speech (Convention Hall)

13:50-14:10

1. Prof. Waonho Yi
INTERNATIONAL COLLABORATIVE NETWORK FOR MT. BAEKDU VOLCANIC DISASTER

14:10-14:30

2. Prof. Hong Huang
STUDY ON VULNERABILITY AND COMPREHENSIVE RISK ANALYSIS METHOD AND SYSTEMDEVELOPMENT BASED ON GIS

14:30-14:50

3. Prof. Chisachi Kato
INDUSTRIAL APPLICATIONS OF LARGE-SCALE FLUID FLOW ANALYSIS

Moderators: Dr. Ho Choi
14:50-15:30

Group photo – Coffee break
Room As311·312

15:30-17:30

Room As303·304

Session 1:

Session 2:

Session 3:

StructuralEngineering (1)Computer Aided Design

Biomechanics andComputer supported
Collaborative Design

Prof. Heecheul Kim
18:00-20:00

Room As313·314

Prof. Jooho Choi

Dr. Yudai Honma

Banquet (Capo Pericano)

CODE 2015

―7―

Time

Tuesday, 15 December
Room As311·312

09:00-10:30

Room As313·314

Session 4:

Session 5:

Session 6:

Disaster Risks

Computational
Mechanics (1)Programing and
Simulation

Dr. Takaaki Kato

Prof. Haengki Lee

10:30-11:00

11:00-12:30

Dr. Muneyoshi Numada

Coffee break
Session 7 :Session 8 :StructuralEngineering (2)Computational
Mechanics (2)Prof. Hong Huang

Dr. Kohei Nagai

12:30-13:45

13:45-15:15

Room As303·304

Lunch
Session 9:

Session 10:

StructuralEngineering (3)Computational
Mechanics (3)Prof. Sangjoon Shin

Prof. Cheekyeong Kim

15:15-15:45

Coffee break
Keynote Speech (Convention Hall)
1. Prof. Sangho Lee

15:45-16:05

DISASTER INFORMATION MANAGEMENT USING 3D CITY INFORMATION MODEL

2. Prof. Marie Oshima
16:05-16:25
16:25-16:45

A NEW CHALLENGE: APPLICATIONS OF SIMULATION TO CLINICAL STUDY

3. Prof. Heoungjae Chun
BIOMECHANICAL ANALYSIS OF SURGICAL SIMULATION FOR DISC DEGENERATION ON THE
ADJACENT SEGMENTS AFTER LUMBAR FUSION AND DECOMPRESSION SURGERY

Moderators: Dr. Hideomi Gokon
Closing Ceremony (Convention Hall)

16:45-17:00

Closing Speech (Prof. Waonho Yi, Dept. of Architectural Engineering, Kwangwoon
Univ., Seoul, Korea)
Closing Speech (Prof. Kimiro Meguro, Director of ICUS, IIS, The University of Tokyo,
Japan)
Moderators: Dr. Hideomi Gokon―8―3rd International Conference on Computational Design in Engineering

Keynote Speech (Convention Hall)

1. Prof. Waonho Yi

INTERNATIONAL COLLABORATIVE NETWORK FOR MT. BAEKDU VOLCANIC DISASTER

2. Prof. Hong Huang

STUDY ON VULNERABILITY AND COMPREHENSIVE RISK ANALYSIS METHOD AND
SYSTEM DEVELOPMENT BASED ON GIS

3. Prof. Chisachi Kato
INDUSTRIAL APPLICATIONS OF LARGE-SCALE FLUID FLOW ANALYSIS - Present status
and future perspectives -4. Prof. Sangho Lee

DISASTER INFORMATION MANAGEMENT USING 3D CITY INFORMATION MODEL

5. Prof. Marie Oshima

A NEW CHALLENGE: APPLICATIONS OF SIMULATION TO CLINICAL STUDY

6. Prof. Heongjae Chun

BIOMECHANICAL ANALYSIS OF SURGICAL SIMULATION FOR DISC DEGENERATION ON
THE ADJACENT SEGMENTS AFTER LUMBAR FUSION AND DECOMPRESSION SURGERY

CODE 2015―9―

TECHNICAL PROGRAMME – DAY 1
PARALLEL SESSION 1:
Structural Engineering (1)Venue: Room As311·312
Date/Time: 14 December, 2015 MONDAY/ 15:30-17:30
15 minutes including 12 minutes presentation and 3 minutes discussion
Session Chair: Prof. Heecheul Kim

1. A. KHAN and H.S. KIM

FREQUENCY RESPONSE ANALYSIS OF SMART COMPOSITE LAMINATES WITH
PARTIALLY DEBONDED PIEZOELECTRIC SENSOR

2. R.NOMURA, S.TAKASE, S. MORIGUCHI and K. TERADA

MULTISCALE ANARYSISWITH A POROSITY MODEL FOR EVALUATING DAMPING
EFFECTSOF COASTALFOREST

3. S. RAJASEKHARAN, M. NUMADA and K. MEGUROSIMPLIFIED NONLINEAR REINFORCED CONCRETE ANALYSIS USING A FINITEELEMENT MAPPED SPRING NETWORK

4. Y. YANG, K. MATSUKAWA, H. CHOI and Y. NAKANOEVALUATION OF RESIDUAL AXIAL LOAD-CARRYING CAPACITY OF RC COLUMNSAFTER SHEAR FAILURE

5. P. JIRADILOK and K. NAGAINUMERICAL SIMULATION OF CHLORIDE DIFFUSION IN CRACKED CONCRETE BASEDON TRUSS NETWORK AND RBSM

6. P.R. DENG and T. MATSUMOTOREBAR FORCE ESTIAMTION FROM CRACK MOUTH OPENING DISPLACMENT

7. C.H. EOM and H.C. KIM

STRENGTH EVALUATION OF HYDRAULIC CEMENT MORTAR WITH GRAPHENE― 10 ―

3rd International Conference on Computational Design in Engineering

PARALLEL SESSION 2:
Computer Aided Design
Venue: Room As313·314
Date/Time: 14 December, 2015 MONDAY/ 15:30-17:30
15 minutes including 12 minutes presentation and 3 minutes discussion
Session Chair: Prof. Jooho Choi

1. J.C. LEE, J.H. JUNG and C.K. KIM
REBAR HUB : 3D REBAR AUTO-PLACING SYSTEM BASED ON BIM

2. J.H. YU and C. KIM
VIBRATION/SOUND
STRUCTURES

REDUCTIONOF

AUTOMOTIVETRANSMISSION

HOUSING

3. H.S. KIM and J.H. CHOIA STUDY ON STATISTICAL VALIDATION OF ELASTO-PLASTIC ANALYSIS OF PISTON
INSERTION INTO HOUSING

4. I. TACHIBANA, S. MORIGUCHI, K. TERADA and S. TAKASEA STUDY ON DIRECT PERMEABILITY SIMULATINON OF GEOMATERIALS USINGIMMERSED BOUNDARY METHOD
5. H.H.S. NAGALUR SUBRAVETI, M. IRYO and T. OGUCHFEASIBILITY STUDY OF INTERACTIVE PASSING MANEUVER MODEL FOR TWO-WAY
TWO-LANE ROADS

6. J.R. CHO, Y.H. KIM and Y.H. YOON  MODELING , HOMOGENIZATION AND FATIGUE LIFE EVALUATION OF FABRIC

BRAIDED REINFORCED RUBBER HOSE

7. J. IMANISHI, K. KASHIYAMA and A. KAGEYAMA

THE 3D MESH MODIFICATION SYSTEM FOR TETRAHEDRON ELEMENTS BASED ON VRTECHNOLOGY

8. M.G. CHO, H.S. SHIN, J.M. CHOI, S.H. YANG and S.M. CHANG
MULTISCALE MODELING AND SIMULATION OF POLYMER NANOCOMPOSITESINCLUDING THE AGGLOMERATED FILLERSCODE 2015― 11 ―

PARALLEL SESSION 3:
Biomechanics and Computer supported Collaborative Design
Venue: Room As303·304
Date/Time: 14 December, 2015 MONDAY/ 15:30-17:30
15 minutes including 12 minutes presentation and 3 minutes discussion
Session Chair: Dr. Yudai Honma

1. M. NUMADA and K. MEGURODISASTER MANAGEMENT SYSTEM “BOSS”- CASE FOR HUMAN RESOURCES SIMULATION –2. H. SONOBE, T. ICHIMURA, P.E.B. QUINAY, K. FUJITA and M. HORI

THREE DIMENSIONAL FINITE ELEMENT GROUND MOTION ANALYSIS USING GPGPU

3. T. SHIMURA, T. KATO, T. EDA and D. KOBAYASHIMULTI-TIMES SIMULATION ANALYSIS OF HIGH RISK DURING EVACUATION IN
A POST-EARTHQUAKE URBAN FIRE

4. R. YABE, Y. HONMA and S. TORIMITUANALYZING EFFECTIVE ROAD BLOCKING FOR IMPROVING EVACUATION TIME FROMTSUNAMI USING TRFFIC FLOW MODEL

5. Y. HONMA, Y. MUNEMASA and K. IMAI
MATHEMATICAL MODEL FOR EQUILIBRIUM AND OPTIMUL LOCATION OF
HOUSINGS AND JOBS WITH CONSTRAINED CAPACITY

6. S. KUBO, R. FUJITA, H. YOSHIDA, M. L. L. WIJERATHNE and M. HORI
EVACUATION SIMULATION CONSIDERING BUILDING DAMAGES

7. G. LING, K. KASHIYAMA and J. MATSUMOTO

A 2D-3D HYBRID MODEL BASED ON STABILIZED FINITE ELEMENT METHOD FORTSUNAMI RUNUP SIMULATION

8. X.Z. LU, X. ZENG, T.Y. YANG and Z. XU

APPLICATION OF THE FEMA P-58 METHOD FOR REGIONAL SEISMIC LOSS
PREDICTION― 12 ―

3rd International Conference on Computational Design in Engineering

TECHNICAL PROGRAMME – DAY 2

PARALLEL SESSION 4:
Disaster Risks
Venue: Room As311·312
Date/Time: 15 December, 2015 TUESDAY/ 9:00-10:30
15 minutes including 12 minutes presentation and 3 minutes discussion
Session Chair: Dr. Takaaki Kato

1. Y. BHATTACHARYA and T. KATODEVELOPMENT OF A VOLUNTARY DECISION-MAKING MODEL FOR RESIDENTS
DURING POST-DISASTER RECOVERY

2. M. OTA, G. LING and K. KASHIYAMASTABILIZED FINITE ELEMENT METHOD BASED ON VOF METHOD FOR FREE SURFACEFLOW USING LES

3. J.P. ZHANG, J.G. CHEN and H.Y. YUAN

PRE-ESTIMATION METHOD FOR LANDSLIDE HAZARD WITH MONITORING ANDSIMULATION

4. T. SRITARAPIPAT and W. TAKEUCHI

PREPARING SAFE URBAN ROADS FOR FLOOD AND EARTHQUAKE IN YANGON,
MYANMAR5. T. KOTANI, S. TAKASE, S.MORIGUCHI, K. TERADA,Y. FUKUTANI, Y.
Otake, K. NOJIMA and M. SAKURABA

PROBABILISTICANALYSIS OF TSUNAMI HAZARDUSINGNUMERICAL SIMULATIONSANDRESPONSESURFACE

6. H. GOKON, S. HAYASHI, S. KOSHIMURA and K. MEGUROTOWARDS AN EXPECTATION OF TSUNAMI FEATURES USING SAR DATACODE 2015― 13 ―

PARALLEL SESSION 5:
Computational Mechanics (1)Venue: Room As313·314
Date/Time: 15 December, 2015 TUESDAY/ 9:00-10:30
15 minutes including 12 minutes presentation and 3 minutes discussion
Session Chair: Prof. Haengki Lee

1. J.W. JUNG, H.W. JANG and J.W. HONG
NUMERICAL SIMULATION OF HIGH VELOCITY IMPACT ON STEEL-PLATE CONCRETE
WALL

2. S. TAKATSU and M. NUMADAINFORMATION SYSTEM "COCOA" FOR EFFECTIVE EVACUATION OPERATION BYPERSONAL IDENTITY NUMBER

3. R.OKAMURA, H. YOSHIKAWA, T. TAKAHASHI and K. KASHIYAMAA ROAD TRAFFIC NOISE SIMULATION METHOD BASED ONACOUSTIC WAVE
THEORYUSINGTIME DOMAIN FM-BEM

4. J.M. CHOI, J.H. MOON and M.G. CHO
A METHOD OF MOLECULAR COMPUTATION OF PHOTO-RESPONSIVE POLYMERREINFORCED WITH NANOPARTICLES

5. S.HORIIKE and K. KASHIYAMASTABILIZED FINITE ELEMENT METHOD FORWIND FLOW ANALYSISWITHPLANT
CANOPY6. J.Q. HAN, J.G. CHEN, H.Y. YUAN and J.P. ZHANGNUMERICAL INVESTIGATION OF NATURAL CONVECTION EFFECTS COUPLED WITH
THERMOACOUSTIC OSCILLATING FLOW IN DIFFERENT GRAVITY INCLINATION― 14 ―

3rd International Conference on Computational Design in Engineering

PARALLEL SESSION 6:
Programing and Simulation
Venue: Room As303·304
Date/Time: 15 December, 2015 TUESDAY/ 9:00-10:30
15 minutes including 12 minutes presentation and 3 minutes discussion
Session Chair: Muneyoshi Numada

1. J. LEE and S.W. YOONSTRUCTURAL TOPOLOGY OPTIMIZATION OF ELECTROMAGNETIC VIBRATION ENERGYHARVESTERS

2. B. SU, H. HUANG and Y. LI

OIL SPILL SPREAD SIMULATION MODEL ON COMPLEX TOPOGRAPHY BASED ON
HYDRODYNAMICS

3. H.Y. CHUNG, J.S. PARK, J.M. CHOI, J.H. YUN and M.G. CHO

A FINITE ELEMENT STUDY ABOUT MECHANICAL FRUSTRATION OF NEMATIC SOLIDS
INDUCED BY DISCLINATION DEFECT

4. M. IRYO and W.K.M. ALHAJYASEENPEDESTRIAN SPEED CHANGE MODEL FOR SAFETY ASSESSMENT AT SIGNALIZEDCROSSWALKS

5. E.Y. KIM and M.G. CHO

NON-INTRUSIVE REDUCTION METHOD WITH CONNECTED ELEMENTS

6. G.H. YOONSTRESS CONSTRAINTS FOR DUCTILE AND BRITTLE MATERIALS IN FLUID-STRUCTURE
INTERACTION SYSTEMCODE 2015― 15 ―

PARALLEL SESSION 7:
Structural Engineering (2)Venue: Room As311·312
Date/Time: 15 December, 2015 TUESDAY/ 11:00-12:30
15 minutes including 12 minutes presentation and 3 minutes discussion
Session Chair: Prof. Hong Huang

1. H. SATO, M. KOSHIHARA and T. MIYAKEOPTIMAL DESIGN OF TIMBER THROUGH COLUMN

2. Y. PATEL and K.MEGURO

NUMERICAL ANALYSIS OF 2D-INFILLED FRAME WITH INTERFACE MATERIALSBETWEEN REINFORCED CONCRETE FRAME AND INFILL WALL

3. B.J. YANG, HAMMAD, R. KHALID and H.K. LEE

INFLUENCES OF ENTRAPPED AIR VOIDS ON THE MECHANICAL AND STRENGTHENING
CHARACTERISTICS OF SPRAYED FIBER-REINFORCED POLYMERIC (SFRP)COMPOSITES: MICROMECHANICS-BASED PARAMETRIC ANALYSIS

4. D.H. SHIN, S.J. HONG and H.J. KIM

EQUATION FOR PREDICTIVE HORIZONTAL RESPONSE SPECTARA IN KOREA

5. M. SAFDAR, T. MATSUMAOTO and K.KAKUMAFINITE ELEMENT ANALYSIS OF REINFORCED CONCRETE BEAMS REPAIRED WITHULTRA-HIGH PERFORMANCE FIBER REINFORCED CONCRETE (UHPFRC)

6. T. ASAI, K. MATSUKAWA, H. CHOI and Y. NAKANORESPONSE ESTIMATION METHOD OF REINFORCEDCONCRETE BUILDINGS DUE TO
WATERBORNE DEBRIS IMPACT LOADS― 16 ―

3rd International Conference on Computational Design in Engineering

PARALLEL SESSION 8:
Computational Mechanics (2)Venue: Room As313·314
Date/Time: 15 December, 2015 TUESDAY/ 11:00-12:30
15 minutes including 12 minutes presentation and 3 minutes discussion
Session Chair: Dr. Kohei Nagai

1. J.W. KANG, H.D. Lee and B. KIM
FULL-WAVEFORM INVERSION FOR MATERIAL PROFILE RECONSTRUCTION IN PMLTRUNCATED DOMAINS

2. S.Y. Tang, X.M. Shu, S.F. Shen and J. Hu

STUDY ON PORTABLE CELLPHONE DETECTOR FOR EARTHQUAKE

3. Z.Y. KANG, T.Q. BUI and S. HIROSECRACK PROPAGATION SIMULATION BY THE
INTERPOLATION FINITE ELEMENT METHOD (XCQ4)

EXTENDEDCONSECUTIVE-

4. L.EDDY and K. NAGAISIMULATION OF FAILURE PROCESS OF REINFORCED CONCRETE BEAM COLUMNJOINT WITH MECHANICAL ANCHORAGES BY 3D RBSM

5. K.W. KIM, M.K. KIM, W.Y. CHO and K.Y. YEAUFINITE ELEMENT ANALYSIS OF BURIED PIPELINE ON FORST REGION CONSIDERINGFROST HEAVE

6. C.K. KIM and S.E. LEE

STRAUTO: PARAMETRIC STRUCTURAL MODELLING AND OPTIMIZATION SYSTEMCODE 2015― 17 ―

PARALLEL SESSION 9:
Structural Engineering (3)Venue: Room As311·312
Date/Time: 15 December, 2015 TUESDAY/ 13:45-15:15
15 minutes including 12 minutes presentation and 3 minutes discussion
Session Chairs: Prof. Sangjoon Shin

1. C.K. GADAGAMMA and K. MEGURO
NUMERICAL EVALUATION OF DAMAGE STATES AND SEISMIC CAPACITY OF
REINFORCED CONCRETE COLUMNS

2. A.A.M. DRAR and T. MATSUMOTO

DEVELOPING A NUMERICAL METHOD FOR FATIGUE LIFE PREDICTION OF RC SLABSREINFORCED WITH PLAIN BARS

3. S. EDALO and T. MATSUMOTO

FATIGUE ANALYSIS OF FROST DAMAGED RC SLABS SUBJECTED TO MOVING LOAD
BASED ON BRIDGING STRESS DEGRADATION CONCEPT4. J.S. KIM, T.H. LEE and T.S. HAN
EFFECT OF THE BASE ISOLATION SYSTEM ON THE LNG STORAGE TANK USINGPROBABILISTIC METHOD

5. T.M. KIM, D.H. KIM, M.K. KIM and Y.M. LIM

MONITORING OF A FULL SCALE PRESTRESSD CONCRETE GIRDER DURING
CONSTRUCTION PROCESS

6. I.TSUWA, T.NAKAGAWA and M.KOSHIHARA

3D STRUCTURAL ANALYSIS MODEL AND EARTHQUAKE RESPSONSE ANALYSIS OFWOODEN FIVE-STORIED PAGODA― 18 ―

3rd International Conference on Computational Design in Engineering

PARALLEL SESSION 10:
Computational Mechanics (3)Venue: Room As313·312
Date/Time: 15 December, 2015 TUESDAY/ 13:45-15:15
15 minutes including 12 minutes presentation and 3 minutes discussion
Session Chairs: Prof. Cheekyeong Kim

1. K. PARK, H. T. KIM and D.J. KIM
GENERALIZED FINITE ELEMENT FORMULTATION ON FIBER BEAM ELEMENTS FOR
DISTRIBUTED PLASTICITY ANALYSIS

2. K. YOON, P.S. LEE and D.N. KIM

FINITE ELEMENT FORMULATION FOR FUNCTIONALLY GRADED BEAM STRUCTURES

3. T. YOSHIMACHI, M. TANIGAWA and K. KASHIYAMAACOUSTIC WAVE SIMULATIONBASED ONPSEUDOIMPULSE RESPONSE

4. K.T. KANG, S.H. LIM and H.J. CHUNEVALUATION OF CONTACT AREA ON KNEE JOINT CARTILAGE USING IN-VIVO
SUBJECT-SPECIFIC FINITE ELEMENT MODEL

5. H. CHO, H.S. JOO, C.W. PARK and S.J. SHINNONLINEAR FLEXIBLE MULTIBODY APPROACHES USING CO-ROTATIONAL FINITEELEMENTS

6. H.C. NOH
APPLICATION OF TRANSFORMED AUTOCORRELATION FOR STRUCTURAL SYSTEMSWITH LOGNORMALLY DISTRIBUTED RANDOM VARIABLESCODE 2015― 19 ―

CODE2015

KEYNOTESINTERNATIONAL COLLABORATIVE NETWORK FOR MT. BAEKDU VOLCANIC DISASTER
Waon-Ho Yi1 and Sungsu Lee2
Professor, Dept. of Architectural Engineering, Kwangwoon Univ., Seoul, Korea,
Professor, School for Civil Engineering, Chungbuk National University, Cheongju, Korea,
Correspond to Prof. Waon-Ho Yi (whyi@kw.ac.kr)
1

2

Keywords: Volcanic Disasters, Online Networks, International Expert Database, International Collaboration, VDPRC
1. VOLCANIC DISASTER PREPAREDNESSRESEARCH CENTER
The volcanic eruption in Iceland in 2010 had put airtransport throughout Europe on hold for almost twomonths. Another series of eruptions of Mt. Shinmoedake,in the island of Kyushu, southern Japan, occurred in
2011. In addition, the possibility of Mt. Baekdu eruption
again has been presented since 2002 so that national
interest in volcanic eruptions has recently been on therise. With all of these, the demands for volcanic disaster
response measures are being discussed.
However, little research has been carried out in Korea
on volcanic disasters so far. As such, there are only ahandful of experts in this field, and the focus of researchmanpower and capacity of related fields are scatteredthroughout Korea that the consolidation of the capacity
to respond to volcanic disasters as a varying complex
catastrophe is pressing.
Accordingly, “Volcanic Disaster Preparedness
Research Center” was launched for R&D projects ofdisaster preparedness technologies since 2012.
These projects of the center is expected to gather
prompt information about volcanic activities, enhance
the capacity to convey the information nationwide,
evaluate national volcanic disaster risks to mitigate
damages from the disasters, and consolidate the ability torespond to them. It is also to solve early on the problem
of a lack of experience and technology to respond to thisnew type of disaster, evaluate its risks, and consolidate
the technologies to respond to them for a safe society.
2.
INTERNATIONALCOLLABORATIVENETWORKS FOR VOLCANIC DISASTERS
The purpose of the project is to vitalize an online
network and an expert group network that specializes involcanic disaster response as an international
collaboration system and help the early introduction ofadvanced technologies. For this, the project is to build aninternational expert group network and organize joint
research to build a network base for international
collaboration.
For online networks, 1) similar networks both inKorea and international communities are investigatedand analyzed, 2) an online forum (VDPedia) is
established to respond to volcanic disasters with themeasures drawn up, 3) the data bank of foreigntechnologies and materials related to volcanic disasters
and a system to convey information, are built.
For building databases, 1) Korean and global
technologies specializing in volcanic disasters areinvestigated, 2) technologies held by related foreign
institutes are investigated and databases are built, 3) a
database of international experts is built.
For international collaboration, 1) a framework for aninternational academic congress on how to respond to
volcanic disasters is prepared, and 2) developing
countries of frequent volcanic disasters was collaborated
with.

Fig. 1 Scope of International Collaborative Network
3. CONCLUSIONSThe project produced an online work field, called
VDPedia which is a website incorporating cyber forum
for experts in volcanic disasters and database showing
information on volcanic damage such as area, casualties,
and etc. Firstly, VDPedia is providing the online space
making experts in volcanic activity, volcanic disaster,and emergency management discuss on the activity andinduced disaster. Its goal is establishing the international
cooperation system on the volcanic disaster and it makes
us quickly response to recover the devastated district. It
is also available to secure technical components via
various expertise establishments and its sharing.
Furthermore, VDPedia is providing current status on
the volcanic disaster with a particular emphasis on Asia.
The information is mainly focusing on volcanic activityinduced damages.
At last, VDPedia hopes to be a hub for continuous
research and discussion on response to the volcanic
disaster.
4. ACKNOWLEDGEMENT
This research was supported by a grant [MPSS-NH2015-81] through the Natural Hazard Mitigation
Research Group funded by Ministry of Public Safety andSecurity of Korean government.

― 21 ―

STUDY ON VULNERABILITY AND COMPREHENSIVE RISK ANALYSIS METHOD AND
SYSTEM DEVELOPMENT BASED ON GIS
Hong HUANG1, Boni SU2, Jinlong ZHAO2, and Nan ZHANG2
1

Professor, Institute of Public Safety Research, Department of Engineering Physics, Tsinghua University, Beijing, China
2
Institute of Public Safety Research, Department of Engineering Physics, Tsinghua University, Beijing, China
Correspond to Prof. Hong HUANG (hhong@tsinghua.edu.cn)

Keywords: Vulnerability, Integrated risk analysis, Multi-hazard, Disaster chain, GIS
1. INTRODUCTION
Increasing various risk sources and factors in cities
make public incidents complicated and diversified.
Multi-hazard and disaster chain are typicalcharacteristics and big challenges. As prevention is thepremise for public safety, the vulnerability analysis andcomprehensive risk analysis are highly required.
In this study, the method of the combination ofstochastic analysis and dynamic simulation weredeveloped to study the vulnerability analysis andcomprehensive risk analysis in urban area, especially
focusing on the multi-hazard and disaster chain. Then, an
integrated system based on GIS was developed which
has been applied in the emergency department of several
cities.
2. METHODS
An emergency management concept based on thetriangle theory of public safety was proposed by Fan [1].
The theory considers that the whole process of theemergency and its response includes three key aspects:
the emergencies, the acceptors, and the emergency
management.
In this study, a comprehensive risk assessment model
was established based on the triangle theory of public
safety. The model considers risk from three aspects:
hazardous of emergencies, vulnerability of acceptors, andcapacity of emergency management. The system takes all
the main factors of public safety into account, and lays agood foundation for the comprehensive risk assessmentof urban public safety.
Recently, multi-hazard and disaster chain become
typical characteristics and big challenges to urban public
safety. The multi-hazard and disaster chain includes
typhoon, heavy rain, water-logging, debris, dangerous
gas leakage, explosion, fire, traffic congestion and so on.
Such emergencies may cause great casualties andeconomic losses. It is important to improve thequantitative risk analysis method of disaster chains.
In this study, the evolution of disaster chain issimulated and the effect’s region and level is analyzed
for risk analysis. Take typhoon for example, the structure
diagram of disaster chain risk assessment model isshown in Figure 1. Typhoon brings rainstorm, rainstorm
cause regional landslide and waterlogging, waterlogging
leads to traffic congestion. The disaster chain model
integrates multiple disaster simulation models, including
typhoon wind field simulation model, typhoon rainfall

simulation model, rainfall-runoff simulation model, andtraffic simulation model. The intensity and character ofdifferent emergency events were simulated using thesemodels. Therefore, spatial-temporal dynamic risk
analysis could be achieved.
Figure 1. Structure diagram of disaster chain risk
assessment model
Landslide
History Data
RainfallDistribution
TyphoonRainfall Model
TyphoonData

DEM Data
Drainage
System Data

RainfallLandslide Model
TyphoonDisaster
Chain Model

Waterlogging
Simulation
Model

Land Use
Data

Traffic
Simulation
Model
Water
Distribution

DEM Data

Geology
Data
Land Slide
Risk

Road Network
Data

Traffic Situation

An integrated system based on GIS was developed toapply the risk analysis method. The main functions of thesystem are comprehensive risk map, dynamic risk
assessment, emergency response capability analysis, andevent risk management. The system has been used inseveral cities and achieved good results.
3. CONCLUSIONSIn this study, the method of vulnerability analysis andcomprehensive risk analysis for multi-hazard anddisaster chain was developed. An integrated system
based on GIS was developed to apply the risk analysis
method to actual cities. The method and system provide
a good approach for public safety assurance.
REFERENCES[1] W. Fan, Y. Liu, Urban Management Science &Technology. 11(2009) 38-41.

― 22 ―INDUSTRIAL APPLICATIONS OF LARGE-SCALE FLUID FLOW ANALYSIS
Chisachi Kato
Professor, Institute of Industrial Science, The University of Tokyo
Meguro-ku, Tokyo, Japan
Keywords: Large-scale simulation, Industrial Application, Computational Fluid Dynamics, Aeroacoustics
For the predictions of turbomachinery flows and flowrelated phenomena that take place in turbomachinery suchas vibrations and noise, the steady and/or unsteady RANS
(Reynolds-Averaged Navier-Stokes) equations have been
widely used in the industries. Since they model the effectsof turbulent eddies in all the scales, simulations based on
RANS have an inevitable limitation for their prediction
accuracy, which has been preventing them from becoming
a complete alternative for testing. Recently, hybrid
methods such as DES (Detached-Eddy Simulation) thatcombines RANS near wall and LES (Large-Eddy
Simulation) elsewhere has been proposed and applied tothe predictions of industrial flows such as those inturbomachinery and those around an automobile. DES
relies on, however, RANS model near the wall andtherefore the prediction accuracy for near-wall velocity
profiles and/or separation from the wall, which is criticallyimportant for turbomachinery flows, does not essentially
go beyond that of RANS.
On the other hand, with tremendous increase inperformance of supercomputers, the range of applicationsof computational methods to predict fluid flows and flowrelated phenomena are spreading rapidly. In fact, as ofNovember 2014, the theoretical peak performance ofhigh-end supercomputers in the world are in several tens
peta-flops and it is expected to reach one exa-flops inaround 2020. Fluid-flow computations that use several
tens billion grids has already become feasible and the
maximum number of the computational grids is expected
to reach a trillion even for industrial applications in 2020.
With that number of the computational grids, fullyresolved LES will become feasible for almost all the typesof turbomachinery flows. Since fully-resolved LES
directly computes smallest but active eddies in the flows,
it can provide as accurate results as DNS (Direct
Numerical Simulation) does. The only difference between
the fully-resolved LES and DNS is that the former models
the effects of eddies that are responsible for dissipation
while the latter computes them directly with thecomputational grids.
Our ultimate goal is to provide industrial users withgeneral purpose simulation software, or to be more exact,
software system capable of performing fully-resolved
LES for industrial flows. To that end, we have beendeveloping FrontFlow/blue system. This presentation will

give a brief introduction to the system, and explain
numerical methods and particular features of the system.
The latest applications that uses several tens billion grids
will also be presented. Finally, our perspectives for thecomputations of the turbomachinery flows for 2020 andafterwards will be given and plans for the near-future
projects will also be explained.

― 23 ―DISASTER INFORMATION MANAGEMENT USING 3D CITY INFORMATION MODEL
S.-H. LEE1 and S.I. PARK2
Professor, School of Civil & Environmental Engineering, Yonsei University, Seoul, Korea,
Ph.D. Candidate, School of Civil & Environmental Engineering, Yonsei University, Seoul, Korea
Correspond to Prof. S.-H. Lee (lee@yonsei.ac.kr)
1

2

Keywords: Disaster Information Management, City Information Modeling, Information Model
1. INTRODUCTION
Development of information modeling technology havebecome an environment in which there are no majordifficulties in expressing 3D objects with attributes. The
approaches should have a different method in case of city
information modeling compared to an individual structure.
This study discusses the considerations to utilize the city
information model for managing the disaster information.
2. CITY INFORMATION MODELING FOR
DISASTER INFORMATION MANAGEMENT
The building and the city information modeling have asimilarity in aspects of 3D geometric model
implementation and management of semantic meanings
by attributes. The strategy for city information modeling,
however, should be different from building information
modeling by following characteristics; the city
information model includes numerous individual models
while BIM focuses the only one structure, the city
information model have to consider the terrain features,and the objectives of the city and building applications aredifferent.
Level of Development (LoD). Modeling is the processthat abstracts required information from the real world.
The required information depends on the tasks ofapplications. The LoD guides the development steps forgeometric features and their corresponding semantic
attributes. The low level of LoD model guarantees theeffective simulation process in case of city information
model-based application which contains huge data, while
the high level model can provide detailed information.
The relationships between facilities. The city facilities
are linked closely. For example, the height of terrain
influences skyline of buildings directly, the roadconditions are the one of factors to estimate the evacuationroute. Therefore, it is indispensable part to consider the
relationship between city facilities. The authors were
proposed the terrain-centered identifiers to link theinformation of facilities according to the features thatevery city facility is located on terrain [1]. Especially,
relationship elements in data schema for city information
modeling can guarantee the consistency of connection
method [2].
Management of dynamic data. City information
model-based applications require additional dynamic datasuch as temperature or water level according to the task ofsimulation. The dynamic data should be managed throughgeneric elements in data schema [2]. Actual dynamic datastored in external database can be connected with the city
information model by generic ones. These features can
support the reuse of the city information model in otherapplications. Figure 1 represents the city information
modeling concept based on the CityGML data schema.
Figure 2 shows an example of flooding simulation using
city information model.
3. CONCLUSIONS
The authors emphasized the LoD of 3D models, therelationship between city facilities, and management of
external dynamic data to build the city information modelfor disaster information management. It is a significant
issue to support the disaster information management
using 3D city information model. Therefore, concrete
application processes based on the city information model
are required through future studies.
Figure 1. City information modeling concept

Figure 2. An example of model-based flooding
simulation (CityGML LoD1 model with hydraulic data)
ACKNOWLEDGEMENT
This work was supported by the National Research
Foundation of Korea(NRF) grant funded by the Korea
government(MSIP) (No. 2011-0030040).
REFERENCES[1] S.I. Park, M.-S. Kim, J.M. Kim, and S.-H. Lee (2015),
City information model-based information managementof flood damages, J. Comput. Struct. Eng. Int. Korea,
28(4), 385-392 (in Korean).
[2] S.-H. Lee, S.I. Park, J. Park, and K.-W. Seo (2014),
Open BIM-based information modeling of railway
bridges and its application concept, ICCCBE 2014, 504511.

― 24 ―

A NEW CHALLENGE: APPLICATIONS OF SIMULATION TO CLINICAL STUDY
M. OSHIMA1, M.KOBAYASHI2, and H. ZHNAG3
Professor, Interfaculty Initiative in Information Studies, The University of Tokyo, Tokyo, Japan,
Ph.D Student, Interfaculty Initiative in Information Studies, The University of Tokyo, Tokyo, Japan,
3
Ph.D Student, Graduate School of Engineering, The University of Tokyo, Tokyo, Japan,
Correspond to Professor M. OSHIMA (marie@iis.u-tokyo.ac.jp)
1

2

Keywords: Multi-Scale Simulation, Medical Image Data, Circulatory System, The Circle of Willis,
Carotid Artery Stenting, Carotid Endarterectomy
1. INTRODUCTION
A technique of simulation has been widely used invarious fields of engineering. Nowadays, it has been
expanding its horizon to the medical and clinical studies.
The authors have been developing an integrated
simulation system combined with medical image data toinvestigate the effects of the surgery such as carotid arterystenting (CAS) or carotid endarterectomy (CEA) on the
cerebrovascular circulation for a patient with a serious
stenosis due to atherosclerosis. Since widening thestenosis area affects not only the diseased area but also theentire blood circulatory system, it is necessary to consider
the entire circulatory system. Diameters of blood vessels
in the circulatory system vary in orders of magnitude frommm to m. In order to simulate the entire circulatory
system, the paper developed the multi-scale numericalmethod based on the one dimension (1D) - zero dimension(0D) [1] or the three dimension (3D) – 1D – 0D model [2].
The present method was applied to the patients withstenosis and the results were compared to themeasurement data.
2. NUMERICAL METHOD
Figure 1 shows the schematic illustration of thepresent multi-scale simulation system.

original 1D-0D simulation for the entire circulatory
system. The governing equations for the 1D model are theequations of mass conservation, momentum conservation,and the relationship between the pressure and the area
while those for the 0D model are the equations of massand momentums conservation. Finally, the results arevisualized in the 3D by remapping the 1D results onto thepatient-specific 3D vascular geometry.
When the patient-specific geometry was applied to the
1D-0D simulation, a deviation from the original 1D-0D
simulation [1] rises. In order to minimize the deviation,
the single photon emission computed tomography(SPECT) data were applied as the reference data to theefferent arteries. Three cases were simulated. The results
showed the good agreement with the phase-contrast (PC)MRA data.
3. CONCLUSIONSThe paper presented an integrated numerical system of
the 1D-0D simulation with the entire circulatory system
for an individual patient using the patient data from the
medical image data. The numerical method wasdeveloped to adjust the differences between the original
1D-0D and the resulting flow rates using the SPECT data.
The present method was applied to three different patientsand the results collated well with the measurement data.
ACKNOWLEDGEMENTS
This study is supported by the Grant-in-Aid for
Scientific Research (B) (No.24300158) by Ministry of
Education, Culture, Sports, Science and Technology inJapan. The authors would like to thank Dr. Shigeli Yamada(Rakuwakai Otowa Hospital) for the medical image data,
and Professor Fuyou Liang (Shanghai Jiao Tong
University, China) and Professor Shu Takagi (The
University of Tokyo, Japan) for their variable discussion
on the 1D-0D simulation.

Figure 1: The schematic illustration of the multi-scale
simulation system
First of all, the patient vascular geometry of the circleof Wills is extracted from the computer tomography (CT)or the magnetic resonance angiography (MRA) data using
our modeling program, V-Modeler [3]. Second of all, thepatient-specific geometric model was incorporated into

REFERENCES[1] F. Y. Liang, S. Takagi, R. Himeno, H. Liu, Medical
Biol. Eng. Comput., 47(2009)743-755.
[2] M. Oshima, R. Torii, S. Tokuda, S. Yamada, and A.
Koizumi, Curr. Pharm. Biot. 13, (2012) 2153-2165.
[3] M. Kobayashi, K. Hoshina, S. Yamamoto, Y. Nemoto,
T. Akai, K. Shigematsu, T. Watanabe, and M.
Ohshima, Circulation Journal, doi: 10.1253/circj.CJ15-0037, 2015.

― 25 ―

BIOMECHANICAL ANALYSIS OF SURGICAL SIMULATION FOR DISC DEGENERATION ON
THE ADJACENT SEGMENTS AFTER LUMBAR FUSION AND DECOMPRESSION SURGERY
H.J. CHUN1 and K.T. KANG1
1Intelligent Structure and Integrated Design Laboratory, Department of Mechanical Engineering, Yonsei University,
Seoul, Republic of Korea
Correspond to Prof. H.J. CHUN (hjchun@yonsei.ac.kr)

Keywords: Disc degeneration, Lumbar fusion surgery, Lumbar spinal stenosis, Spinous process osteotomies,
Conventional laminectomy, Adjacent segment disease, Finite element analysis
1. INTRODUCTION
The presence of degeneration alters the stress withinthat particular degenerated disc [1]. The state ofdegeneration of the disc adjacent to a single level fusion
and after lumbar decompression surgery would influence
the rate of subsequent further degeneration of that sameadjacent disc. Therefore, the purpose of this study was toinvestigate the biomechanical effect of the discdegeneration at the adjacent segment on the stress of thecorresponding segment after one segment lumbar fusion
and two different lumbar decompression surgeries using
validated Finite Element (FE) models.

percent change at L3-L4 from intact to decompression
models was less in the moderate disc degeneration modelcompared to healthy disc model (Fig 2 and 3).

2. MATERIAL AND METHOD
Three validated L2–L5 FE models with different
grades of disc degeneration (normal, mild, moderate) at
the L2–L5 were simulated. Based on the intact models
that represented each grade of disc degeneration, thematched
L3–L4
instrumented
fusionanddecompression using conventional laminectomy (CL)and spinous process osteotomy (SPiO) models weredeveloped. Loading and boundary conditions werereferred to our previous study using spine FE model [2].

Figure. 1. The effect of disc degeneration and fusion on
maximal von Mises stress of AF and intradiscal pressureof NP at L2–L3 and L4–L5 segments under hybrid
testing protocol

3. RESULT
One-level lumbar fusion surgery
Disc degeneration led to an increase of maximal vonMises stress of annulus fibrosus at this degenerated disc,
related to disc degeneration under flexion, extension,lateral bending and torsion moments. Furthermore, thefusion also increased a maximal von Mises stress ofannulus fibrous at both proximal and distal adjacentsegments under hybrid protocol. The increases of annular
stress were greater at the proximal adjacent segment thanthose at distal adjacent segment. However, thecontribution of the fusion surgery to the increment rate ofannulus fibrosus stress at the adjacent segment did notincrease in relation with disc degeneration states at theadjacent segment (Fig. 1).
Conventional laminectomy and spinous process
osteotomy lumbar decompression surgery
The CL model and moderate disc degeneration model
demonstrated a larger annulus stress under all four
moments than those in the SPiO model and healthy discmodel, respectively. Therefore, the decompression level(L3-L4) in the CL model with moderate discdegeneration showed the highest annulus stress. The

Figure 2. The maximal von Mises stress in the AF of theintact, (a) SiPO and (b) CL models related to discdegeneration under four moments

Figure 3. NP in the intact, (a) SiPO and (b) CL models
related to disc degeneration
4. CONCLUSIONSIn conclusion, disc degeneration per se and lumbar
fusion and decompression can cause the increase of discstress at the adjacent segment. However, they seem not
to act synergistically. In addition, although the CL modelin the moderate disc degeneration showed the highest
annulus stress, the degenerative models would be lessinfluenced by the decompression technique.
REFERENCES[1] L. Bresnahan, A.T. Ogden, R.N. Natarajan, R.G.
Fessler, Spine. 34 (1976) 17–23.
[2] H.J. Kim, H.J. Chun, K.T. Kang, H.M. Lee, B.S.
Chang, C.K. Lee, J.S. Yeom, Yonsei Med J. 56 (2015)
146-153― 26 ―

CODE2015

Technical Program – Day 1

Session 1

FREQUENCY RESPONSE ANALYSIS OF SMART COMPOSITE LAMINATES WITH
PARTIALLY DEBONDED PIEZOELECTRIC SENSOR
A. KHAN1, and H.S. KIM2
Research Assistant, Dept. of Mechanical, Robotics and Energy Eng., Dongguk University, Seoul, Republic of Korea,
Associate Professor, Dept. of Mechanical, Robotics and Energy Eng., Dongguk University, Seoul, Republic of Korea,
Correspond to Prof. H.S.KIM (heungsoo@dgu.edu)

1
2

Keywords: composite laminates, piezoelectric sensor, debonding, layerwise theory, frequency response
1. INTRODUCTION
In the last few decades, composite laminates bonded orembedded with smart materials, named as smart
composite laminates, have been developed and widelyused in space vehicles, fixed-wing aircraft, rotary-wing
aircraft, civil structures, marine systems, automotive
systems, robotics systems, machine tools, and medical
systems [1]. Piezoelectric material are the most popularsmart materials. In smart composite laminates,
piezoelectric sensors and actuators are used for active
vibration
control,delamination
detectionandclassification, structural health monitoring, etc.
Literature review reveals that in most of theapplications, the piezoelectric patches are considered to be
perfectly bonded to the host structures [2]. However, thedebonding failures between the host structures andpiezoelectric elements may occur during the service life,
due to bonding defects, impact loading or high peelingstress concentration at the bonding edge. Debonded sensor
may results in an erroneous sensing signal which may
adversely affect the proposed applications. Therefore, it isimportant to evaluate the effects of partial sensor
debonding on the dynamic behavior of the smart structures.
2. MATHEMATICAL FORMULATION
Piezoelectric sensors and actuators are considered asadditional layers of the composite laminate. Improved
layerwise theory with Heaviside unit step function [3] isemployed for the displacement field, while the potential
variation through the thickness of piezoelectric patches ismodeled using higher order electric potential field. The
coupled electro-mechanical field with discontinuities atthe debonded interface, is implemented with finite
element method. The governing equation is obtained
through the Extended Hamilton’s principle. The
developed model is numerically implemented on a smart
composite plate with one piezoelectric actuator and one
piezoelectric sensor.
3. FREQUENCY RESPONSE ANALYSIS
To obtain the frequency response of the smart
composite plate with partially debonded sensor, the plateis excited by an electric harmonic load with the magnitudeof 1V. General mode superposition method is applied tosolve the 2nd-order governing differential equation in
frequency domain. For partial debonding at the right sideedge of the sensor, the decrement in the magnitude of eachpeak in the frequency response curves, normalized withrespect to perfectly bonded sensor case, is shown in the

Figure 1. It is observed that the magnitude of sensor output
at the resonant frequency has reduced almost linearly withincreasing length of debonding.

1
1st natural frequency
2nd natural frequency
3rd natural frequency

0.95
0.9
0.85
0.8
0.75
0.7
0.65
0.6

0

5

10

15

20
25
30
Debonding length (%)

35

40

45

50

Figure 1. Effect of partial edge debonding in the sensor on
the magnitude of peaks in the FRF of sensor output
3. CONCLUSIONS
Partial sensor debonding at the edge of piezoelectric
sensor in smart composite laminates results in thedecrement of the magnitude of sensor output at theresonant frequencies. The shift in resonant frequencies is
subtle.
ACKNOWLEDGEMENT
This work was supported by the Technology Innovation
Program (10048305, Launching Plug-in Digital Analysis
Framework for Modular system Design) funded by the
Ministry of Trade, Industry & Energy (MI, Korea).
REFERENCES[1] I. Chopra., J. Sirohi, Smart Structures Theory,
Cambridge University Press, 2013.
[2]H. Boukabache, C. Escriba, J.Y. Fourniols, Toward
smart aerospace structures: Design of a piezoelectric
sensor and its analog interface for flaw detection. Sensors
2014, 14, 20543-20561.
[3]H.S. Kim, A. Chattopadhyay, A. Ghoshal,
Characterization of delamination effect on composite
laminates using a new generalized layerwise approach.
Computers & structures 2003, 81, 1555-1566.

― 27 ―

MULTISCALE ANARYSIS WITH A POROSITY MODELFOR EVALUATING DAMPING EFFECTS OF COASTAL FOREST
R. NOMURA1, S. TAKASE2, S. MORIGUCHI3 and K. TERADA4
1
Student, Graduate school of Engineering, Tohoku University, Sendai, Japan,
Assistant Professor, Department of Civil and Environmental Engineering, Tohoku University, Sendai, Japan
3
Associate Professor, International Research Institute of Disaster Science, Tohoku University, Sendai, Japan
4
Professor, International Research Institute of Disaster Science, Tohoku University, Sendai, Japan
Keywords: Coastal forest, Multiscale analysis, Tsunami run-up, Porous media, Stabilized finite element method2

1. INTORODUCTION
A coastal forest is often built for disaster reduction
against tsunami. The quantitative evaluation of its
damping effect is required to predict the inundation area
and the damages on structures behind the forest.
Although several studies have been conducted on thissubject of study, most of them are based on the hydraulicexperiments. Therefore, it is still difficult to estimate thedamping effect quantitatively in real situations. To solve
this problem, we propose a multiscale numericalapproach that enables us to evaluate the damping effectof the coastal forest in a wide area with relatively low
computational costs. The proposed method also allows
us to evaluate the influence of the fluid force acting on
structures. In this study, a simple numerical example isconducted to demonstrate the performance of the
proposed method.
2. NUMERICAL ANARYSISThe proposed multiscale numerical analysis method is asort of homogenization methods and designed tocharacterize the damping effect of coastal forests byperforming numerical analyses of the meso- and microstructures. More specifically, a coastal forest is stratified
into three spatial scales as shown in Figure 1; the branchscale (micro), the tree-scale (meso) and the forest-scale(macro) and the lower-scale analysis is expected toprovide the damping effect for the upper-scale analysis.
To transmit the lower-scale information about thedamping effect to the upper-scale, a porosity model is
employed. That is, by regarding the coastal forest as aporous media in the forest-scale, we expect that theoverall damping effect can be evaluated by themechanical behavior at tree-scale.
In numerical analyses, fluid flow is simulated with thestabilized Finite Element Method (FEM) developed tosolve the Navier-Stokes equation accompanied by aporosity-related drag term [1], which characterizes thedamping effect.

Figure 1. Multiscale analysis: three spatial scales

3. COMPUTATIONAL EXAMPLE
Figure 2 shows the result of the numerical experimentat the forest-scale. Two calculation cases are considered;
Cases with (a) low porosity and (b) high porosity. Eachcase has a porous domain that represent the coastal forest(Length=160m, Width=140m). A building is also set
behind the coastal forest. Based on results of the branchscale and the tree-scale analyses, a coefficient of theporous term is defined as a function of the fluid velocity
and depth. As we can see from Figure 2, the fluid forCase (b) reaches 180m longer than that of Case (a). This
result indicates that a dense forest is superior to a sparseone in resistance against fluid flow and consistent withthe outline of the result obtained from hydraulic experiment. Figure 3 shows the result of the stress distribution
on the surface of the building obtained from structural
analysis.

Figure 3. Result of numerical experiment (time=60s)

Figure 4. Stress on the building (time=6.3s)
4. CONCLUSION
We propose the framework of multiscale analyses forevaluating the disaster-reduction effect of coastal forests.
The presented multiscale numerical method enables us toevaluate not only the damping effect but also the stresses
in structures caused by the fluid force.
REFERENCESG[1] F. Aprino, N. Massarotti, A. Mauro. International
Journal of Heat and Mass Transfer 54 (2011) 4495–4506.

― 28 ―

SIMPLIFIED NONLINEAR REINFORCED CONCRETE ANALYSISUSING A FINITE ELEMENT MAPPED SPRING NETWORK
S. RAJASEKHARAN1, M. NUMADA2, and K. MEGURO3
Project Researcher, Institute of Industrial Science(IIS), The University of Tokyo(UTokyo), Tokyo, Japan,
2
Lecturer, IIS, UTokyo, Tokyo, Japan,
3
Professor and Director, International Center for Urban Safety Engineering(ICUS), IIS, UTokyo, Japan
Correspond to S. RAJASEKHARAN (shan@iis.u-tokyo.ac.jp)

1

Keywords: reinforced concrete, smeared crack model, spring network, secant formulation
1. INTRODUTION
There are many uncertainties involved in the seismic
vulnerability assessment of a large existing building stock,which involves a large number of simulations, to obtain areliable probabilistic vulnerability estimate. In developing
countries, the presence of a large number of nonengineered reinforced concrete (RC) buildings, further
adds to the level of uncertainty involved and thecomputation effort required.
To this effect, a simplified method to analyze the nonlinear behavior of RC is proposed in this research. Thismethod involves the usage of a finite element mappedspring network [1], using spatially averaged material
models [2]. A secant stiffness based formulation is used[3]. The RC domain is discretized into an assemblage ofconcrete and steel springs.
2. EXPERIMENTAL VALIDATION
The applicability of the model is checked bycomparing the numerical simulation with experimentallytested RC panels [4], RC beams [5], (Figure 1) and RC
frames [6] (Figure 2). Good correlation was obtained withthe experimental data of the RC beams and RC Frames,
using only the concrete compression/tension model. The
applicability of concrete shear models is verified throughthe experimental validation of RC panels subjected toshear loading. Better shear models are required in
specimens where excessive shear slippage occours.

Through the spring network the crack pattern in an
averaged sense can also be obtained.

(a) Experimental
Setup(b) Spring
Network(c) Observed
Crack
Pattern

Figure 2. Experimental validation of a RC frame
3. CONCLUSIONS
The results show the capability of the proposed model
to follow the nonlinear behavior of simple RC structuralmembers while using a reduced number of degrees of
freedom, which also reduces the computation effort
required. Simplified spring stiffness formulation anddiscretization are used, which can be efficiently used tofollow the nonlinear behavior of simple structuralelements. The advantage of the proposed spring networkis that it can be used in tandem with a connected discreteelement model, which can be used to perform futher large
deformation analysis.

(a) Experimental Setup(b) Spring Network(c) Observed Crack Pattern

REFERENCES[1] A.A. Gusev, Phys. Rev. Lett. 93, 034302 (2004).
[2] H. Okamura, K. Maekawa, Nonlinear analysis andconstitutive models of reinforced concrete, Gihodo Co.
Ltd., Tokyo, 1991.
[3] F. Vecchio, Amer. Concrete Inst. Strut. J., 86 (1989),
26-35.
[4] F. J. Vecchio, M. P. Collins, J. Amer. Concrete Inst.,
83-22 (1986), 219-231.
[5] B. Bresler, A.C. Scordelis, J. Amer. Concrete Inst., 60(1963), 51-72[6] K. Muto, Strength and deformation of structures,
Maruzen, Tokyo, 1965.

Figure 1. Experimental validation of a RC beam― 29 ―EVALUATION OF RESIDUAL AXIAL LOAD-CARRYING CAPACITY OF RC COLUMNSAFTER SHEAR FAILURE
Y. YANG1, K. MATSUKAWA2, H. CHOI2, and Y. NAKANO3
Graduate Student, Department of Architecture, The University of Tokyo, Tokyo, Japan,
Research Associate, Institute of Industrial Science, The University of Tokyo, Tokyo, Japan,
3
Professor, Institute of Industrial Science, The University of Tokyo, Tokyo, Japan,
Correspond to Y. YANG (yangyong@iis.u-tokyo.ac.jp)
1

2

Keywords: Shear Failure, Arch Resistance Model, Residual Axial Load Carrying Capacity

3. APPLICATION OF THE ARCH RESISTANCE
MODEL
To verify the accuracy of the evaluation formula ofresidual axial load-carrying capacity, a database of test
results is compiled from the previous studies. TheN

Q=0

NS QS MSNC
QC
arch effect

δ

e

NC
QC

NS QS MS

Figure 1 Model of shear-damaged RC columns(1)nN S   2nM S  N C e
2000

1000

0

0
1000
2000
Measured axial load (kN)(a) Arch resistance modelCalculated axial load (kN)

2. ARCH RESISTANT MODEL
In this paper, the state of shear force equal to 0 isdefined as the limit state of axial collapse. The model ofshear-damaged RC columns can be established as Figure
1. The equilibrium equation of moment of it is shown inEquation 1.
It should be noted that in this model the axial forces
acting with an eccentricity e at the end sections ofconfined concrete part can develop a force couple, which
resists the moment induced by the P-Δ effect of
longitudinal steel bars. The resistance of the force couple(Nce) acting on the crushed concrete can be considered as
the interaction between crushed concrete and longitudinal
steel bars and it is called ‘arch effect’ in this research. It isthe big difference compared with other models proposedin the past studies.
Based on the theory of structural mechanics, theevaluation formula of the residual axial load-carrying
capacity of shear-damaged RC columns can be deduced.

application results of the shear-friction model [2] and thearch resistant model with the same database are shown in
Figure 2. It shows that the arch resistant model can give a
better estimate than the shear-friction model for most of
specimens.

Calculated axial load (kN)

1. INTRODUCTION
From the post-earthquake reconnaissance survey, it isobserved that reinforced concrete (RC) short columns andRC columns with poor transverse reinforcement arevulnerable to shear failure. For severely shear-damaged
RC columns without adjacent load redistribution members
around them, the deterioration of the axial load carrying
capacity of these damaged columns can lead to a partial or
global building collapse. It can cause a great loss of lives
and property. Therefore, it is necessary to evaluate theresidual axial load-carrying capacity of them.
Until now, several evaluation models [1]-[3] have been
proposed by other researchers to account for the loss of
axial load-carrying capacity for such RC columns prone
to shear failure. However, it is not easy to obtain an
intuitive understanding regarding the loss of axial load
carrying capacity through these models.
Hence, in this paper, an arch resistance model is
proposed for shear damaged RC columns, which can givea better understanding of the loss of axial load-carrying
capacity.

2000

1000

0

0
1000
2000
Measured axial load (kN)(b) Shear-friction model

Figure 2 Comparison of calculated-to-measured results
4. CONCLUSIONSFor RC columns damaged in shear, the arch resistancemodel to predict residual axial load-carrying capacity is
proposed. For most of the specimens included in thecompiled database, the calculated residual axial loadcarrying capacity through the arch resistance model has agood agreement with the measured result. Its highaccuracy is confirmed by comparing with the application
result of the shear-friction model with the same database.
REFERENCES[1] Y. Uchida, Y. Uezono, Journal of Structural andConstruction Engineering. 566 (2003) 177-184. (inJapanese)[2] K.J. EIwood, J.P. Moehle, ACI Structural Journal. 102(2005) 578-587.
[3] Y. Takaine, M. Yoshimura, Journal of Structural andConstruction Engineering. 618 (2007) 191-197. (inJapanese)― 30 ―

NUMERICAL SIMULATION OF CHLORIDE DIFFUSION IN CRACKED CONCRETE BASEDON TRUSS NETWORK AND RBSM
Punyawut JIRADILOK 1 and Kohei NAGAI 2
1

Doctoral course student, Department of Civil Engineering, The University of Tokyo, Japan
2

Associate Professor, Institute of Industrial Science, The University of Tokyo, Japan

Keywords: chloride ion penetration, truss network model, RBSM, crack width, Fick’s 2nd law
1. GENERAL AND BACKGROUND
Many reinforced concrete structures are exposed to themarine environments which result in chloride ion
penetration; subsequently decreasing in structural
capacity due to the corrosion. Therefore, an estimation ofchloride ion penetration in concrete with cracking isimportant for maintenance and predicting service life of
the structures. This study aims to develop a threedimensional numerical simulation system of chloride ion
penetration in cracked concrete based on truss network
and Rigid Body Spring Model (RBSM).
2. ANALYSIS MODEL
In this study, the truss network system is developed andincorporated into RBSM to simulate the chloride ion
penetration in sound and cracked concrete. Truss networkmodel is developed based on Fick’s 2nd law for simulating
the time dependent chloride ion transportation, while
RBSM can simulate cracks from the interaction of eachconcrete element at meso-scale which calculated by the
respond of spring that set on interface between
elements.[1] In practice, the chloride ion transport throughcrack concrete is much faster than that of sound concrete,because of additional flow channel provided by crack
space. In RBSM simulation, when concrete is loaded, thecrack propagation is directly induced by the spring
between elements, and chloride diffusion coefficientthrough the crack spacing increase according to the crack
width which equal to spring elongation [2].
3. GEOMETRY OF NUMERICAL SIMULATION
Chloride exposed

Line load

(b)

(a)

Fixed

Chloride exposed surface
Crack geometry in 3D

Figure 1. (a) Boundary condition for RBSM numericalloading, (b) Chloride exposed boundary condition.
The indirect tensile loading test model is conducted inorder to obtain the uniform crack shape. The shortcylinder numerical model size of 100 mm diameter and
50mm thickness are subject to indirect tensile loading. A
line of compressive displacement is set at one side ofnumerical model, while the opposite side is set a line of
fixed movement condition. The load was given until thedesignated crack width passed through specimens (40 ,80 , 100 ). Then, the chloride exposure condition is
given at one flat surface of cracked concrete models for
simulate the chloride penetration in cracked concrete.
3. RESULTS

Experiment

Chloride profile by simulation

Figure 2. Comparing chloride contamination profile withexperiment.
Chloride penetration profiles from simulations werecompared with previous experiment by K. Audenaert(2009)[3] and show a good agreement. The chloride
transportation through wider crack width is much fasterthan that of narrow crack width and sound concrete part,the chloride concentration at the perpendicular zone to thecrack direction is also higher.
4. CONCLUSION
The coupling simulation can correctly predict thechloride penetration profile and represent the effect ofdifferent crack geometry appropriately comparing withthe experiment studies.
REFERENCES[1] K.Nagai, Y. Sato and T.Ueda. Mesoscopic simulationof failure of mortar and concrete by 3D RBSM.Journal ofAdvanced Concrete Technology, 2005, pp.385-402.
[2] A. Djerbi, S. Bonnet, A. Khelidj and V. BaroghelBouny, Influence of traversing crack on chloride diffusion
into concrete , Cement Concrete Res, 2008, pp.877-883.
[3] K. Audenaert, G. De Schutter and L. Marsavina V.
Boel, Influence of cracks on chloride penetration andcorrosion initiation time, 2nd International RILEM
Workshop on Concrete Durability and Service Life
Planning, 2009.

― 31 ―

REBAR FORCE ESTIAMTION FROM CRACK MOUTH OPENING DISPLACMENT
P.R. DENG1 and T. MATSUMOTO2
1

Ph.D. Student, Graduate School of Engineering, Hokkaido University, Sapporo, Japan,
2
Associate Professor, Faculty of Engineering, Hokkaido University, Sapporo, Japan

Keywords: Crack mouth opening displacement (CMOD); Fracture mechanics; Reinforced concrete; Bond slip
1. INTORDUCTION
For reinforced concrete (RC) structures, the state ofrebar is an important indicator of structural performance.
Fortunately, the rebar force has an inherent relation withconcrete cracking process. This relation has beensuccessfully exploited in rebar force estimation by inverse
analysis of crack opening displacements (COD profile) in[1]. However, COD profile is difficult to be measured inreal existing structures. Therefore, innovation of a reliable
and convenient technique deserves attention.
In this paper, the transformation to relate CMOD andrebar forces of RC structures under applied loads isestablished based on fracture mechanics and bond slip
model. Correspondingly, a rebar force estimation method
is proposed by using the easily measurable data, CMOD.
2. PROBLEM FORMULATION
According to the concrete crack mechanics, the total
CMOD (um) of an RC beam is divided into three parts:
CMOD due to applied loads (uma), rebar forces (umb) andbond slip (ums). The relation is expressed as:um  um a  umb  um s
(1)among which uma and umb can be calculated based on
fracture mechanics and are given by:
4 a  a(2)
um a 
G x , a , b  x dx G x, a , b da 
E  0  0 

4 a  a
(3)
G x , a , b  f x dx G x, a , b da 
E  0  0
where E’ is a combination of elastic constants; a is thecrack length and a’ is the dummy variable of a; G isweight function which can be found in [2]; σ(x) and f(x)are stress that would exist on the crack face in the absenceof a crack due to applied load and rebar forces per unit
length, respectively.
The ums can be calculated by employing a local bondslip model in [3]. Consequently, the ums is expressed witha function of rebar force f(x) as: a
 f x dx um s  S    S  0
(4)
 db





where ε and db are rebar strain and diameter, respectively.
um a   

3. APPLICATIONIn this paper, experimental COD data of a tested beam
under different load conditions in [1] is employed inestimating the rebar force of that beam. The experimental
CMODs are listed in Table 1.

Table 1. Experimental CMOD
Load (kN) 17.64 19.60
CMOD(mm) 0.132 0.156

21.56
0.192

23.52
0.221

25.48
0.257

Table 2. Experimental CMOD
Load
(kN)
17.64
19.60
21.56
23.52
25.48

Rebar force (kN/mm2)
(1)

(2)

(3)

304.2
257.8
318.3
290.4
314.5

244.6
278.4
264.1
278.0
287.3

217.9
242.1
266.2
280.4
318.3

Estimating error(%)[(3)[(3)(1)]/(3)(2)]/(3)-39.59-12.23-6.51-15.00-19.56
0.44
-3.57
4.26
1.48
8.66

Table 2. shows the rebar force estimation and thecorresponding errors from different approaches: (1)Inverse analysis of COD profile by employing method in[1]; (2) CMOD analysis method proposed in this paper; (3)Standard RC cracked beam transformed section analysis,
the results of which can be treated as theoretical rebar
force because the bridging effect of concrete is negligible
comparing with that of rebar force and the test data is froma newly cured beam under four-point static bending test.
It is found in Table 2 that the rebar force estimating
accuracy of CMOD analysis method is generally higher
than the accepted method, inverse analysis of COD profile.
4. CONCLUSION
Following the mechanisms of concrete cracking for RC
structures, a rebar estimation method from CMOD basedon fracture mechanics and bond slip model is proposed inthis study. The proposed method is successfully employedin estimating rebar force of the cracked beam in [1]. The
better accuracy than the generally accepted rebar force
estimation method, inverse analysis of COD profile,verifies the applicability and reliability of this method.
REFERENCES[1] I. Nazmul, T. Matsumoto, International Journal of
Solids and Structures. 45 (1) (2008) 159-174.
[2] H. Tada, P. Paris, G. Irwin, The Analysis of Cracks
Handbook, ASME press, New York, 2000.
[3] K. Maekawa, H. Okamura, A. Pimanmas, Non-linear
Mechanics of Reinforced Concrete, CRC press, 2003.

― 32 ―

STRENGTH EVALUATION OF HYDRAULIC CEMENT MORTAR WITH GRAPHENE
Chul-Hwan EOM1, Heecheul KIM2
1

Professor, Department of Architecture & Fire Protection, Uiduk University, Gyeongju, South Korea,
2
Professor, Department of Architectural Engineering, Kyung Hee University, Yongin, South Korea,
Correspond to Prof. Heecheul KIM (kimhc@khu.ac.kr)

Keywords: Carbon Nano Tube Composite material, Graphene, Mortar, Strength enhancement
1. INTRODUCTION
Graphene’s extraordinary properties such as about 200
times stronger than steel by weight shows possibility of
excellent role in enhancing the strength of a concrete.
Graphene composite was introduced as a strength
enhancement material referencing the research on carbon
nano tube composite since the material properties of bothmaterials are very similar. Strength improvement effect ofgraphene composite was observed both on compressive
and tensile fields.
2. EXPERIMENTAL PLAN
Strength evaluation experiment was performed
according to the ASTM regulation. Mixing ratios ofwater-cement and nano materials which would affect boththe compressive and tensile strength were consideredduring mixing process of mortar. The quantity of nanomaterial was also considered as one of the importantparameters during mixing process.
Table 1. Characteristics of Graphene
Test Items

Specification

Analysis
Results

Units

C-content
Ash
Moisture
Bulk density

>99.5
<0.5
<0.5
/

99.64
0.36
0.23
145

%
%
%
g/l

3. EXPERIMENT RESULT AND ANALYSIS
The water-cement ratio being kept constant through all
specimen. The quantity of graphene was selected as amain parameter for the strength change measurement. The
quantity of graphene was 0 to 20% of weighted ratio to thecement.
Figure 1. Compressive
strength experiment &specimen
About 20% compressive
strength increment wasobserved along with the 7%graphene composite mortar as shown in table 1. Also thesimilar percentage of tensile strength was observed alongwith the 5~7% graphene composite mortar.

Figure 2. Tensile strength experiment & specimen

Table 1. Results of test
Strength(MPa)Compressive
Tensile
GM00
36.6
2.0
GM01
42.4
2.4
GM02
43.6
2.2
GM03
43.8
1.9
GM04
43.3
2.0
GM05
42.4
2.4
GM06
44.2
2.4
GM07
41.6
2.3
GM08
34.4
1.8
GM09
32.3
2.1
GM10
30.3
1.9
GM15
29.7
2.2
1)
GM20
29.8
1.2
GM201) : 20 : Amount of Graphene (%)
Specimen

4. CONCLUSIONSGraphene contributes to the improvement of bothcompressive and tensile strength of mortar. The ratio of
improvement was 20% in maximum based on the amountof graphene included in mortar.
REFERENCES[1] ASTM C109M-13e1 Standard Test Method forCompressive Strength of Hydraulic Cement Mortars(using 2-in. or [50-mm] Cube Specimens),[2] G.M. Kim, S.C. Choi, C.S. Shim, Y.S. Chung,Compressive strength evaluation of composite mortar
with grapheme-oxide, Proceedings of the KCI Fall
Conference of 2012(in Korean).

― 33 ―

CODE2015

Technical Program – Day 1

Session 2

REBAR HUB : 3D REBAR AUTO-PLACING SYSTEM BASED ON BIM
J.C. LEE1, J.H. JUNG2, and C.K. KIM3
Professor, Dept. of Architectural Eng., Tongmyong University, Busan-si, Korea,
Professor, School of Architecture, Kyungnam University, Changwon-si, Korea,
3
Professor, Dept. of Architectural Eng., Dankook University, Yongin-si, Korea,
Correspond to Prof. C.K.KIM (cheekim@dankook.ac.kr)

1

2

Keywords: BIM, Rebar Detailing, Auto-placing, Structural design
1. INTRODUCTION
A 3-dimensional rebar auto-placing system, Rebar Hub,has been designed and implemented in this research.
Rebar Hub provides an integrated BIM platform from 3D
structural modeling of buildings to detailed rebar autoplacing which reflects precise anchorage lengths, lap
splices, and order lengths. Rebar Hub executes 3D rebar
auto-placing by recognizing 2D drawing CAD files which
are auto-converted into 3D structural models. After the
auto-placing, each structural member composing the 3D
structural model will hold rebar properties data of its own.
It means that these rebar properties data can be used forthe subsequent works such as quantity-survey,
manufacturing, and rebar fabrication.
2. REBAR AUTO-PLACING
Based on the generated 3D frame model, Rebar Hub
auto-generates 3D rebar shapes and placing drawings(Figure 3). Rebar Hub auto-places bars in precise 3D
shapes by referring to the relevant rebar detailing codesuch as KBC or ACI. Precise lap splices and seismic
detailing are reflected in the auto-placing, and othervarious rebar detail options can be set by the user. Unlike
other BIM tools where rebar placing is done one element
by one element, which can be said to be semi-automaticand time-consuming and all the individual structuralmembers have to be handled by hand, Rebar Hub executes
auto-placing for a project by one button-click as all theneeded data have been built from 2D recognition. This is
the key feature that differentiates Rebar Hub from all theother current BIM tools.
From the generated 3D rebar model, Rebar Hub autogenerates 2D shop and placing drawings required for onsite works (Figure 4). 2D rebar detail drawing works havebeen the most time-consuming process in the conventional
project delivery, but this has been removed in Rebar Hub.
Rebar Hub also provides 2D manual editing tools forfrequent on-site modifications.
3. CONCLUSION
The development of Rebar Hub started from this verybottom-up approach. Using Rebar Hub, the work can startfrom conventional 2D drawings, which are converted to
3D, while being able to fully edit 2D figures with the helpof familiar 2D authorizing tools. The application of Rebar
Hub to a number of pilot projects has proved to be
successful so far, to be a feasible BIM tool integrating 2Dand 3D modeling.

Figure 1. 3D rebar auto-placing by Rebar Hub

Figure 2. Placing Drawing Auto-generation

Figure 3. Quantity estimation generated by Rebar Hub
REFERENCES[1] Rundell, R. (2008). "BIM and Digital Fabrication (12-3 Revit Tutorial)."― 35 ―

VIBRATION/SOUND REDUCTION OF AUTOMOTIVE TRANSMISSION HOUSINGSTRUCTURES
J.H. YU1 and C. KIM2
1

Graduate Student, Department of Mechanical Engineering, Kyungpook National University, Daegu, Korea,
2
Professor, Department of Mechanical Engineering, Kyungpook National University, Daegu, Korea,
Correspond to Prof. C. KIM (kimchul@knu.ac.kr)

Keywords: Manual Transmission, Vibration Transfer, Transmission Housing, Bracket Noise, BEM
1. FEM/BEM MODEL
The housing structure of automotive transmission for apassenger car was modeled with many NASTRAN shell
finite elements as in Fig. 1, to analyze forced vibration
responses caused by fluctuations of an internal
combustion engine. The noise is then transmitted fromthe T/M housing through TGS cables to car interior.
2. ACOUSTIC POWER ANALYSIS
Fig. 2 shows the acoustic field surrounding thetransmission housing structure which would be
investigated to calculate noise power levels in the carinterior. The housing was again modeled with theboundary elements using LMS/Virtual Lab. The forced
response results from FE analyses are input to BEM
mesh nodes of housing, and the sound pressure at eachnode is, then, calculated. Finally, SPL can be obtained at
the wanted field points on the acoustic field surrounding
the housing model. It is known from Fig. 2 that the soundpressure levels are high at the back (driver position) of
the housing, up to 136 dB in the car interior.
The acoustic power levels felt from a driver and apassenger positions are shown in Fig. 3. The upper line
denotes the acoustic power level at a driver position.
Peak dB level was noticed in the frequency of 2300 Hz.
The difference was caused by the distance between adriver and a passenger.

Figure 2. The acoustic field point mesh and results

3. CONCLUSIONS
Design changes in shapes of brackets and thethickness of supporting rubber washers installed on an
automotive transmission housing, that are based on FEM
and BEM sound and vibration analyses, can effectively
reduce the level of sound and vibration transferred to thevehicle interior.

Figure 3. The acoustic power levels in a car interior
REFERENCES
Reference to a journal publication:[1] K.S. Oh, S.K. Lee, S.J. Kim, KSNVE, 2007(3) (2007)
861-866.
Reference to a book:[2] H.I. Lee, J. Sound and Vibration, 20(7) (2013) 585596.
[3] Y.H. Kown, C. Kim, D.W. Kim, and S.H. Jung,
KSAE, 2008(5) (2008) 639-644.

(a) FE Model of T/M(b) Forced response
Figure 1. A FE model of T/M housing and responses― 36 ―

A STUDY ON STATISTICAL VALIDATION OF ELASTO-PLASTIC ANALYSIS OF PISTON
INSERTION INTO HOUSING
H.S. KIM1, J.H. CHOI2
1

Department of Aerospace and Mechanical Engineering, Korea Aerospace University, Republic of Korea
2

School of Aerospace and Mechanical Engineering, Korea Aerospace University, Republic of Korea,
Correspond to Prof. J.H. CHOI (jhchoi@kau.ac.kr)

Keywords: Statistical Calibration, Validation, Bayesian Approach, Elasto-Plastic Analysis
1. INTRODUCTION.
A Pyrotechnically Actuated Device (PAD) is acomponent that delivers high power in remoteenvironments by the combustion of a self-containedenergy source. Once the PAD triggers combustion in achamber, the expanding gas drives a piston into housing.
It usually plays safety-critical role of valve that opens a
gas flow in aviation and defense industry such as
launching and deployment. Thus the reliability analysis of
the PAD is very important.
In this study, statistical calibration and validation
procedure is addressed for the elasto-plastic insertion
analysis of the piston into the housing by the generated
gas pressure. Objective of the study is to develop reliable
analysis that agrees closely with the actual test data, andvalidate the accuracy in probabilistic manner.
2. MODELING & CALIBRATION.
Toward this objective, the actual test data are provided
for three types of PAD. The result of test is the resistanceforce against the piston movement, which is measured bymaterials testing system under quasi-static condition. In
order to evaluate the resistance force via FEA, commercial
analysis tool ANSYS is used and elasto-plastic solution is
obtained between the piston and housing during the quasistatic piston movement. In the analysis, there are a couple
of issues in view of uncertainty. First, the friction
coefficient is not clearly identified, which can vary withmany factors such as surface condition and contact forces.

whole device. Then the boundary condition of the model
is not clearly specified, but is imposed arbitrarily. All
these are the sources of disagreement of the analysis resultwith the test data. In this case, common practice is tointroduce the associated parameters as the unknowns, andcarry out calibration to improve the accuracy of theanalysis model which is to match the analysis result withthe actual test data. In this procedure, however, if theparameters are calibrated without considering
uncertainties due to the limited test data, errors of test andmodeling, the result can mislead to wrong conclusions.
To solve this problem, a Bayesian approach is proposed
to calibrate the parameters of computation model of PAD
in a probabilistic way. Markov Chain Monte Carlotechnique is employed, which is an efficientcomputational algorithm to draw samples that satisfies theposterior distribution.
3. VALIDATION & CONCLUSIONS.
Once the calibration is carried out for a baseline model,
one obtains the parameters in the form of distribution. The
calibrated input parameters are applied to the analysis ofother types of PAD. The result of analysis is validated inprobabilistic way at each type of PAD. As a means to theprobabilistic validation, area metric is introduced between
the CDF by actual test data and analysis. The moreagreement we get as a result of calibration, the smallerarea we get between the two CDF's. Finally, the analysis
results enclose the test data completely, which means thatthe calibration is made successfully.
REFERENCES[1] Dawn An and Joo-Ho Choi, Improved MCMC method
for parameter estimation based on marginal probability
density function, Mechanical Science and Technology,
Vol . 27, PP. 1771-1779[2] G. Gary Wang, Adaptive response surface method
using inherited latin hypercube design point, American
Society of Mechanical Engineers, Vol. 125, PP.210220

Figure 1. Computational model of piston and housing
Next, the analysis model includes only a part of thedevice as shown in Figure 1, which is the cutout of the― 37 ―

A STUDY ON DIRECT PERMEABILITY SIMULATINON OF GEOMATERIALSUSING IMMERSED BOUNDARY METHOD
I. TACHIBANA1, S. MORIGUCHI2, K. TERADA3, and S. TAKASE4
1

Graduated Student, School of Engineering, Tohoku University, Sendai, Japan,
2
Associate professor, IRIDeS, Tohoku University, Sendai, Japan,
3
Professor, IRIDeS, Tohoku University, Sendai, Japan,
4
Assistant professor, School of Engineering, Tohoku University, Sendai, Japan
Correspond to I. TACHIBANA (tachibana@mm.civil.tohoku.ac.jp)

Keywords: Permeability, Geomaterials, Immersed Boundary Method
1. INTRODUCTION
It still has been difficult to predict detailed response ofgeomaterials theoretically in dynamic interaction
problems between geomaterials and pore fluids, such asliquefaction, erosion, and seepage failure of river dikes.
In practice, on the other hand, there are certain needs ofnumerical simulations focused on geohazards with fluidsolid interaction modeling. In this study, we applied
Immersed Boundary Method (IBM)[1] to a direct
numerical simulation (DNS) on the water passing
through saturated geomaterials. Numerical permeability
tests were conducted and the effects of soil properties,
such as grain size distribution, are discussed.
2. IMMERSED BOUNDARY METHOD
Numerical framework used in this study is composedof Discrete Element Method (DEM), Finite Difference
Method (FDM), and IBM. Soil grains are assumed to be
rigid spheres and the motion of grains is solved by DEM.
On the other hand, flow of pore water is expressed by
FDM using fixed calculation grids. Meanwhile, IBM
accomplishes a conjunction of those two methods.
Here is the Navier-Stokes equation for incompressive
Newtonian fluid with an external force term ݂;డ௨೔డ௧ ൅ ‫ݑ‬௝

డ௨೔

డ௫ೕ

ൌെ

ଵ డ௣

ఘ డ௫೔

൅ߥ

డ

൬

డ௨೔

డ௫ೕ డ௫ೕ൰ ൅ ݂௜ ,(1)

where ‫ ݑ‬is the velocity vector, ‫ ݌‬is the pressure, ߥ isthe kinematic viscosity coefficient, and ɏ is the density.
IBM works to evaluate the external force, which emerges
from the difference of dynamics on the boundaries
between fluid and solid domains, using a smoothed
approximation to the Dirac delta function ߜ௛ [1] forquantity translation between fixed grid points andarbitrary moving points. After the original IBM method
had been proposed, some modifications of thecalculation algorisms have been proposed. In this study,
the one proposed by Uhlmann[2] is used to describeinteraction between soil grains and pore water.
3. PERMEABILITY TESTS
Numerical permeability tests were conducted using 3Dmodel shown in Fig.1. Soil grains (spheres) are fixed inspace. One of the simulated results are shown in Fig.2.
Each plots are with critical Reynolds number [3],indicating there be turbulent-like effects in permeability

flow when it excesses 1. The plot tells us that original
test cases (triangular plots) are getting off the theoretical
proportion in laminar flows as the critical Reynolds
number increases, while they still keep in line when we
modified their pressure conditions to make Reynolds
number lower than 1 for study (circle-plotted cases). It
can be explained supposing there are turbulent-like
effects to diminish flow potential in the soil block.
4. CONCLUSIONS
According to the obtained results, proposed numericalmethod seems to be enough capable to catchcharacteristics of micro-scale flow through granularmaterials. As we hope to describe precise characteristics
of geomaterials in future works, further investigationmust be done in order to have advanced applications.

Fig.1 Soil block substituted by uniformed spheres (left)
and the sliced part of it (right). Pressure gradient ݅ isevaluated with given pressures ܲon in/out planes.

Fig.2 Results of permeability tests with critical Reynolds
number proposed by Fancher, et.al.[3]REFERENCES[1] C. S. Peskin, J. Comput. Phys. 10-2 (1972) 252-271.
[2] M. Uhlmann, J. Comput. Phys. 209 (2005) 448-476.
[3] G. W. Fancher, J.A. Lewis, J.B. Barnes in: Bul. 12,
Mineral In. Exp. St., Penn. State Col., 1933, pp.373-384.

― 38 ―

FEASIBILITY STUDY OF INTERACTIVE PASSING MANEUVER MODEL FOR TWO-WAY
TWO-LANE ROADS
H.H.S NAGALUR SUBRAVETI1, MIHO IRYO2, TAKASHI OGUCHI3
Masters Student, Department of Civil Engineering, University of Tokyo, Tokyo, Japan
Associate Professor, Institute of Industrial Science, University of Tokyo, Tokyo, Japan
3
Professor, Institute of Industrial Science, University of Tokyo, Tokyo, Japan
Correspond to Mr. Nagalur Subraveti (harihs@iis.u-tokyo.ac.jp)

1
2

Keywords: Driving Simulator, Conflict Zone, Action Pattern Choice, Vehicle-to-Vehicle Interaction
1. INTRODUCTION
Vehicular movements on two way two lane roads
differ significantly from other facilities due to the highfrequency of conflicts between the vehicles in opposing
directions of travel. Presence of obstacles like parked
vehicles, construction sites etc. increase the frequency of
conflicts between the vehicles. These conflicts have asignificant impact on the efficiency [1] as well as thesafety of road sections.

Fig 2: Flow chart representing the action pattern process
Fig 1: Representation of the lane with depiction of
conflict zone
Driving Simulators (DS) are used for studying various
aspects of transportation including driving behavior, roadsafety, road design and traffic flow. It is very important tomodel the behavior of surrounding vehicles in the DS
realistically. The existing driver behavior model [2]assumes that the opposing vehicles always have priority
over the passing vehicles. However in reality, drivers of
the passing traffic may adopt aggressive behavior to
obtain priority and the opposing vehicles may giveway
considering such behavior to avoid any potential conflicts.
As it is not possible to deterministically assume thebehavior of the subject driver of the DS, the surrounding
vehicle model should not consider a fixed priority rule but
interactively choose its action. The objective of this study
is to propose a framework of an interactive model between
the passing and opposing vehicles and evaluate thecharacteristics of the model by sensitivity analysis using aMonte-Carlo simulation.
2. ACTION PATTERN
Fig 1 represents a scenario where vehicles travellingin opposing directions interact due to the presence of an
obstacle. Conflict zone is the area where there is apotential for both cars to collide if both assume that thepriority of passing lies with them. Each vehicle has the
choice of either 1) giving way to the other vehicle bydecelerating, 2) employing aggressive movement by
forcing the other vehicle to decelerate to complete themaneuver or 3) move freely without any interaction when

there is a sufficient gap between one’s exit time of theconflict zone and the other’s entering time. The vehicles
are modeled as intelligent agents with the individual goalof reaching its destination as soon as possible and having
the ability to interact and cooperate with other agents tosolve conflicting goals. The flow chart in Fig 2 represents
the action pattern choice to be chosen by the driver. The
vehicles initially check for the feasibility of freemovement and then adopt the aggressive movement incase a collision is expected with the free movement. In thecase of the aggressive movement, the vehicle makes anassumption of the maximum speed decrease the othervehicle is willing to use and check for the feasibility ofcompleting the action with these assumed values. If thisassumption is feasible, the vehicle will take the aggressive
movement by trying to force the other to decelerate,otherwise it gives way.
3. CONCLUSIONSThe proposed model is implemented to a simulation
program. The results of the sensitivity analysis regarding
the impact of model parameters upon the traffic
performance will be shown in the presentation.
4. REFERENCES[1] Durbin, C.T. Traffic performance on two-lane, twoway highways: examination of new analytical approaches.
Diss. Montana State University-Bozeman, College ofEngineering, 2006.
[2] Llorca, C, et al. "Development of a new microscopic
passing maneuver model for two-lane rural roads." Transp.
Res C: 52 (2015): 157-172.

― 39 ―  MODELING , HOMOGENIZATION AND FATIGUE LIFE EVALUATION OF FABRIC

BRAIDED REINFORCED RUBBER HOSE
J.R. CHO1, Y.H. KIM2, and Y.H. YOON3

1

Associate Professor, Department of Naval Architecture and Ocean Engineering, Hongik University, Sejong, Korea,
2
Graduate Student, School of Mechanical Engineering, Pusan National University, Korea,
3
Senior Researcher, Research and Development Institute, MIDAS IT, Gyeonggi-Do,Korea,
Correspond to Prof. J.R. Cho (jrcho@hongik.ac.kr)

Keywords: Fabric braided hose,   modeling, Homogenization, Large deformation analysis, Fatigue evaluation
1. FABRIC BRAIDED BRAKING HOSE INCYCLIC MOTION
Automotive braking hose is in lamination constructionof rubber and fabric braided layers for both oil sealing
and preventing the excessive radial expansion. In
addition, it exhibits the large deformation during thecyclic motion due to the combination of the steering
motion and the up and down motion of vehicle.

Figure 1. Fabric braided braking hose in cyclic motion
The most important function of braking hose is thesafe convey of interior working oil without leakage up to
the predefined mileage of vehicle. However, its reliable
design reliability is not easy because its prediction at thedesign state requires several considerations, such as the -modeling of complex fabric braided layers inserted
between rubber layers, an accurate material modeling forthe reliable large deformation analysis along the specificcyclic path, and the elaborate fatigue model considering
the mean strain.
2.
HOMOGENIZATIONANDFATIGUEEVALUATION
The mechanical behavior of fabric braided layers is
modeled as orthotropic material, and its equivalentmaterial properties were calculated making use of thesuperposition method and the unit cell finite element
analysis [1]. To compute the fatigue life, three fatigue
models, standard model without considering the meanstress and strain, Morrow considering the mean strain
and SWAT(Smith-Watson-Topper) considering the meanstress are employed [2].

Figure 2. Unit cell finite element analysis and the fatigue
life calculation module built in midas NFX.

Meanwhile, the large deformation along the specificcyclic path shown in Figure 1 was carried out by thequadratic interpolation of cyclic path using three known
position and rotation data and by updated Lagrangian
formulation.
3. FATIGUE LIFE PROFILES ALONG THE
SPECIFIC PATH
Figure 3 compares the fatigue life profiles predictedby the proposed homogenization and fatigue evaluation
methods. The left profile is by the standard fatigue model
without considering the initial hose deformation, while
the center and right show the results by Morrow andSWAT(Smith-Watson-Topper) fatigue models. It has
found from the detailed numerical data that theremarkable difference in the critical fatigue lives isobserved even though the overall fatigue life profiles arealmost similar.

Figure 3. Fatigue life profiles of automotive braking hose
for the horizontal cyclic path.
4. CONCLUSIONS
An effective but reliable fatigue life prediction method
was proposed by making use of the homogenization
method and the modified fatigue evaluation models.
5. ACKNOWLEDGEMENT
The research was supported by Basic Science
Research Program through the National Research
Foundation of Korea (NRF) funded by the Ministry of
Education (Grant No. NRF-2014R1A1A2055820).
REFERENCES[1] J.R. Cho, J.I. Song, J.H. Choi, Key Eng. Mater. 306308 (2006) 799-804.
[2] J.R. Cho, Y.H. Kin, C.W. Seo, Y.G. Kim, Finite Elem.
Anal. Des. 100 (2015) 65-76.
[3] J.R. Cho, Y.B. Jee, W.J. Kim, S.R. Han, S.B. Lee.
Composites Part B 53 (2013) 112-120.

― 40 ―

THE 3D MESH MODIFICATION SYSTEM FOR TETRAHEDRON ELEMENTS
BASED ON VR TECHNOLOGY
J. IMANISHI1, K. KASHIYAMA2, and A. KAGEYAMA3
1

Graduate student, Department of Civil and Environmental Engineering, Chuo University, Tokyo, Japan
2
Professor, Department of Civil and Environmental Engineering, Chuo University, Tokyo, Japan
3
Professor, Department of System Informatics, Kobe University, Kobe, Japan,
Correspond to Prof. K.KASHIYAMA (kaz@civil.chuo-u.ac.jp)

Keywords: mesh modification, node relocation, mesh refinement
1. INTRODUCTIONS
The 3D numerical simulations are becoming more
powerful and popular tool for various CAE(Computational Aided Engineering) problems. The
quality of numerical solution depends on the quality ofmesh. If the quality of mesh is extreme worth, we are notoften able to obtain the accurate numerical solution. In
those cases, it is required to modify the mesh. However, it
is difficult to modify the 3D mesh using 2D display.
Therefore, the present authors has been developed amesh modification system for 3D unstructured mesh using
VR technology [1]. Two modification methods are
implemented in this system: the node relocation method
and the mesh refinement method. The previous system has
been shown the efficiency of the method. In order toimprove the applicability of the system, we implement
several new functions such as the node merge function andthe mesh-connectivity modification function.
2. SYSTEM SUMMARY
In order to develop the mesh modification system, we
employed the immersive virtual reality device such as
CAVE. This device is composed by three 3D projectors,
three large screens, position tracking device and active
shutter glasses. The Figure 1 shows the scene that user
modify the mesh in CAVE (VR space) using the system.
User can check the details of 3D mesh and can modify theshape of mesh in VR space interactively by using thecontroller. This system can be applied to the linear
tetrahedron elements and the 2-nd order tetrahedron
elements.
Figure 2 shows the flowchart of this system. After theloading the mesh data, the element quality is computed.
Then, the value of mesh quality is displayed on the screen.
User can perform the mesh modification using theinformation of mesh quality. This system implements two
modification methods: node relocation method and meshrefinement method. For mesh refinment method, thissystem has two refinment patterns (Figure 4). For the noderelocation method, the node control function is
implemented to avoid the distruction of the boundary
surface by the node relocation. This function control thederection of the node movement. Therefore, the element
is modifyed without the destruction of the boundary
surface. We also implements node merge function for theoverlapping elements and mesh-conectivity modification
function for the elements which exist on boundary surface.

Figure 1. VR environments

Figure 2. Flow chart of the system

Figure 3.
Applied elements

Figure 4.
Mesh refinement patterns

The update of mesh data is performed in realtime. Several
application examples will be shown in the presentation.
3. CONCLUSIONSThe 3D unstructured mesh modification system based on
VR Technology for tetrahedron elements has beenpresented. In order to improve the applicability, we have
implemented new functions such as the node merge
function and the mesh-conectivity modification function.
The development of mesh modification for thehexahedron elements is left for the future work.
REFERENCES[1] Tanaka, S., Kashiyama, K., Imanishi, J. andKurumatani, M. Studies on implementation of meshrefinement method and improvement of applicability of
3D unstructured mesh modification system, Journal ofJapan Society of Civil Engineers, Ser. A2 (Applied
Mechanics (AM)), Vol.70 (2014) No.2 pp. I_339-I_348.
(in Japanese)― 41 ―

MULTISCALE MODELING AND SIMULATION OF POLYMER NANOCOMPOSITESINCLUDING THE AGGLOMERATED FILLERS
MAENGHYO CHO1, HYUNSEONG SHIN2, JOONMYUNG CHOI2, SEUNGHWA YANG3, and SEONGMINCHANG4
Professor, School of Mechanical and Aerospace Engineering, Seoul National University, Seoul, South Korea
Ph. D. Student, School of Mechanical and Aerospace Engineering, Seoul National University, Seoul, South Korea
3
Professor, School of Mechanical Engineering, Chung-Ang University, Seoul, South Korea
4
Postdoc researcher, School of Mechanical and Aerospace Engineering, Seoul National University, Seoul, South Korea
1

2

Correspond to Prof. MAENGHYO CHO (mhcho@snu.ac.kr)
Keywords: Agglomeration, Interphase percolation, Molecular dynamics simulation, Finite-element homogenization,
Polymer nanocomposites
1. INTRODUCTION
Filler agglomeration is the most challenging issue in the
manufacturing field of polymer nanocomposites [1]. Dueto the high nonbond interaction such as van der Waals
interaction and electrostatic potential, fillers areagglomerated in the polymeric phase. In spite of theimportance of filler agglomeration, there is no multiscale
framework. In this paper, we propose the multiscalemodeling strategy of polymer nanocomposites including
the agglomerated fillers. To characterize interphasepercolation effect, the molecular dynamics simulations
and finite element homogenization models are employed.
We verify that the proposed interphase percolation model
is applicable to the multiparticulate system.
2. COMPUTATION MODELS AND RESULTS
Figure 1 shows the overview of multiscale modelingstrategy to characterize interphase percolation model andverify the propose model [2]. Spherical silicon carbide
nanoparticles are embedded as the reinforced particle in
the polymer phase (polypropylene). In order to identify
the interphase percolation effect, six molecular models areused for various interphase particulate distance. Choi et al.
[3] proposed the schematics for characterization ofinterphase, which is employed in this study.
As shown in Figure 1, elastic properties of interphasedecreases as the overlapping density increases. The
definition of overlapping density is volume of overlappedzone over the volume of interphase.
Therefore,
overlapped density means how much interphase zone isoverlapped. This is key parameter of this study and it
reflects degradation mechanism of interphase elastic
properties. Interphase percolation model is defined asfollowing form:
overlap
Eint  e

int
Voverlap
int
Vtotal ( overlap   )

1

(1)

 Emat(2)

To verify the extensibility to the multi-particulate
system, additional study is conducted. For same location,
volume fraction, and radii of nanoparticles, thehomogenized elastic properties are compared between

molecular dynamics model and finite element
homogenization models. It shows similar results between
molecular dynamics model (Young’s modulus: 2.55 GPa,and shear modulus: 0.98 GPa), and finite element
homogenization model (Young’s modulus: 2.83 GPa, andshear modulus: 1.07 GPa).

Figure 1. Overview of multiscale modeling strategy toreflect interphase percolation effect
3. CONCLUSIONSInterphase percolation model is proposed and it isdefined by the molecular dynamics models and finite
element homogenization models. Extensibility to themulti-particulate system is verified. The proposed model
will be extended to the elasto-plastic model and toughness
model.
This work was supported by the National Research
Foundation of Korea (NRF) grant funded by the Korea
government (MSIP) (No. 2012R1A3A2048841)
REFERENCES[1] L. Vaisman, H. D. Wagner, G. Marom, Adv. Colloid
Interface Sci. 128 (2006) 37-36.
[2] H. Shin, S. Yang, J. Choi, S. Chang, M. Cho, Chem.
Phys. Lett. 635 (2015) 80-85.
[3] J. Choi, H. Shin, S. Yang, M. Cho, Compos. Struct. 119(2015) 365-376.

― 42 ―

CODE2015

Technical Program – Day 1

Session 3

DISASTER MANAGEMENT SYSTEM “BOSS”- CASE FOR HUMAN RESOURCES SIMULATION Muneyoshi NUMADA1, Kimiro MEGURO2,
Lecturer, ICUS, Institute of Industrial Science, the University of Tokyo, Japan,
Director/Professor, ICUS, Institute of Industrial Science, the University of Tokyo, Japan,
Correspond to Muneyoshi NUMADA (numa@iis.u-tokyo.ac.jp)
1

2

Keywords: disaster response, workload, time chart, schedule, web
1. INTRODUCTION
This research defined the standard response process forthe local governments and developed IT disastermanagement system based on the standard disaster
response process. Developed IT system is named BOSS
(Business Operation Support System or BOSai System,
Bosai means the disaster management in Japanese). BOSS
can show the total disaster process with database ofhistorical disaster knowledge. About 500 processes aredefined by the disaster management plan and theexperiences of staffs in Yabuki town in Fukushima
prefecture and Ishinomaki city in Miyagi prefecture where
was heavily damaged by the 2011 Great East Japan
Earthquake disaster. The core database of BOSS include
people, house, industry, law, land etc. which are defined
the relationship to the disaster responses. Before disaster,
BOSS can be used to simulate disaster responses
according to the level of “scenario disaster”. After disaster,
BOSS can manage the actual responses and work as aplatform of disaster information sharing.
2. Disaster management system “BOSS”
The damage conditions such as number of total collapse
houses, number of death, number of injury, etc. are input
data set of BOSS. The condition of weather, season, time
and week day are also considered.
For each disaster response, the workload is calculated byfollowing to the level of damages. The total amount ofwastes by damaged houses, possible total evacuee,necessary toilets for the evacuee and possible temporary
houses etc. are also calculated. These calculation isderived from the direct damage such as the number ofcollapsed houses.
Human resources are assigned to each disaster response
following each calculated workload. When the humanresources are assigned, the disaster responses are
categorize four types. (1) the responses by everyone canwork immediately after disasters such as management of
the Volunteer Center etc., (2) the responses by government
staff at first stage of disasters but gradually shifting to theother people such as food supply, water supply etc. (3) the
responses by every government staff such as residents
support, issuance of certificate, etc., (4) the responses bygovernment staff with a special skills such as restorationof infrastructure and lifelines, health issues etc.
It is important to make a category for each disaster
response according to these kinds of types. Then we canunderstand what kinds of disaster responses are necessary

by staffs with special skills or without those, and manage
the limited resources effectively. In this paper, to assign
the human resources, three cases can be simulated for
different human resources level (Figure 1).

Figure 1: Comparison of different human resources fordisaster responses
3. CONCLUSIONSBOSS can be used the simulation of disaster responses toshow the critical path and bottle neck response according
to different cases of input conditions. BOSS can show (1)Evaluation of amount of work (work load), (2) Effective
distribution of human resources with skill level and work
load and (3) Management of response schedule. Then
stakeholders can understand all disaster processes withpriority, difficulty or level for the effective humanresource management.
Small local governments cannot prepare enough budget
and knowledge to develop and manage their original IT
system. BOSS is usable for local government to prepare
the disaster response simulation in Japan.

― 43 ―

THREE DIMENSIONAL FINITE ELEMENT GROUND MOTION ANALYSIS USING GPGPU
H. SONOBE1, T. ICHIMURA2, P.E.B. QUINAY3, K. FUJITA4, and M. HORI5
Student, BS, Department of Civil Engineering, University of Tokyo, Japan.
Dr. Eng., Assoc. Prof., Earthquake Research Institute, University of Tokyo, Japan.
3
Dr. Eng., Assis. Prof., Research Institute for Natural Hazards and Disaster Recovery, Niigata University, Japan.
4
Dr. Eng., Visiting Researcher, Advanced Institute for Computational Science, RIKEN, Japan.
5
Ph.D., Prof., Earthquake Research Institute, University of Tokyo, Japan.
Correspond to H. SONOBE (h-sonobe@eri.u-tokyo.ac.jp)
1

2

Keywords: Ground Motion Analysis, Conjugate Gradient Method, GPGPU
1. INTRODUCTION
Large-magnitude earthquakes can severely damage urban
areas, and in order to prepare for seismic hazard of urban
structures, a reliable ground motion analysis technique isessential. To gain accuracy, several aspects need to beconsidered carefully. It is straightforward to consider thephysics of these problems and solve the linear or nonlinear
solid wave equations. Also, the target domain should belarge enough to neglect the boundary effects. Last,
structural damage mainly depends on a frequency band of
up to 2.5 Hz, and the surface topography is known toaffect the response of structures significantly. All of thiscombined leads to a tremendous amount of computationalcost, the degrees of freedom (DOFs) in the order of 100
million to 1 billion, so several research has been
conducted to develop sophisticated computational
technique. One approach is to use CPU based
supercomputers, such as the K computer[1]. Anotherapproach is to use super computers with accelerators, themost popular being GPGPU, general-purpose GPU.
Applying GPU acceleration to seismic analysis is recently
considered to lead superior performance compared to
CPU based computers. [2] uses the finite-different method(FDM), which is capable of gaining high performance dueto the structured data sets. When it comes to modeling thesurface topography accurately, however, FDM needsmaller element sizes, so as a result the DOFs get higher.
[3] uses the spectral-element method (SEM) which uses
higher ordered elements compared to finite-element
method (FEM). In this case, a FE solver would be moresuitable, since it can greatly represent the surface
topography using lower ordered elements, which is easierto formulize and rather simple to implement. Thus, wehave developed a three dimensional finite element groundmotion calculation technique using a heterogeneous
CPU/GPU platform.
2. METHODOLOGY
For this research, we have enhanced the performance ofthe pre-conditioner of a finite element based seismic wave
amplification CG (conjugate gradient) solver, GAMERA,
which was originally developed for CPU based
supercomputers. GAMERA uses a CG solver inside thepre-conditioner to calculate the preconditioning matrix ofthe main CG loop (Figure 1). We have developed aheterogenous computing method for the pre-conditioner,

using linear/second ordered tetrahedral mesh and single
precision arithmetic. We have conducted numericalexperiments to show good weak scaling efficiency, andproved that it is considerably faster compared to using a
CPU cluster of the same size. These data will be shown in
the presentation. The code was rewritten using OpenACC,
a programming standard for parallel computing, and is
known to be a relatively user-friendly way to implementon scientific calculation. The rewritten code was merely
5~10% of the original code, showing high portability as
well as good performance.

Figure 1. Overview of the CG solver in GAMERA and the
proposed method.
3. CONCLUSIONSThis method will be able to enhance forward analysisof seismic response, and treat problems that have higher
DOFs using more computational resources. It will also be
fairly easy to apply to many-core computing architecture,
which is expected to be the mainstream of the next
generation’s supercomputers.
REFERENCES[1] T. Ichimura, K. Fujita, S. Tanaka, M. Hori, M. L. L.
Wijerathne, Y. Shizawa and H. Kobayashi, Physics-based
urban earthquake simulation enhanced by 10.7 BlnDOF
x30 K time-step unstructured FE non-linear seismic wave
simulation, SC '14, 2014.
[2] D. Michéa and D. Komatitsch, Accelerating a threedimensional finite-difference wave propagation codeusing GPU graphics cards, Geophys. J. Int., 2010.
[3] D. Komatitsch, D. Göddeke, G. Erlebacher, D. Michéa,
Modeling the propagation of elastic waves using spectral
elements on a cluster of 192 GPUs, Comput. Sci. Res.
Dev., 2010

― 44 ―

MULTI-TIMES SIMULATION ANALYSIS OF
HIGH RISK DURING EVACUATION IN A POST-EARTHQUAKE URBAN FIRE
T. SHIMURA1, T. KATO2, T. EDA3, and D. KOBAYASHI4
Project Researcher, Institute of Industrial Science, University of Tokyo, Tokyo, Japan,
Associate Professor, Institute of Industrial Science, University of Tokyo, Tokyo, Japan,
3
Researcher, Institute of Industrial Science, University of Tokyo, Tokyo, Japan,
4
Research Student, Institute of Industrial Science, University of Tokyo, Tokyo, Japan,
Correspond to Assoc. Prof. T. KATO (kato-t@iis.u-tokyo.ac.jp)

1
2

Keywords: Disaster Risks, Urban Fire, Evacuation, Agent-based Simulation
1. INTRODUCTION
After the Great East Japan Earthquake in 2011,
initiatives have focused on eliminating the possibility ofan unforeseen disaster situation. In regards to the dangers
of urban fire caused by earthquake in Tokyo, according to
the earthquake damage assumption by the Tokyo
Metropolitan Government in 2012, fatalities are estimatedto about 4,100 people in a fire. It can be said that a highrisk has been estimated.
However, this estimation was performed with aregression analysis using data of past earthquake-induced
fires and peacetime conflagrations all of which occurredin different periods and under different conditions. Thismethod needs to be verified to entail a more scientificapproach. It also does not reveal a real image of thedamage in the current city which has changed
considerably over the years. Therefore, it is necessary toreview again the risk of deaths from urban fire due to an
earthquake in current Tokyo. Additionally, there is also aneed for a better understanding of possible circumstances
for mass mortality like in the case of Great Kanto
Earthquake in 1923.
Therefore, we attempted to understand whether there ismuch possibility of mass mortality by urban fire in themodern metropolis. Furthermore, we consider where thecause of the phenomenon lies and how to avoid it[1].
2. METHODOLOGY
By logic, since people walk at faster pace than the rateof fire spread, simple calculation results often estimatethat people are able to evacuate in time. However, thesecalculations do not take into account that there arepossibilities of crowd congestions which can slow down
the pace of the evacuees thereby entrapping them in fire.
Hence in this research, we developed a simulator whichruns the fire-spread simulation and the evacuation
simulation simultaneously over a wide area[2].
By performing a large number of trials, we haveattempted to clarify the possibilities of mortality due toextreme events and delayed evacuations thereby revealing
the structure of the phenomenon.
Accordingly, the calculation target area wasconcentrated in this range, with the evacuation
destinations considered were particular locations that aredesignated by the Tokyo Metropolitan Government intheir wards.

Hamada type sequence formula was used for the
calculations of the fire simulation model, due to the factthat it employs a high fire-spread rate and it is considered
to be suitable for the evaluation of safety. Occurrence
locations of the fire was elected by a random number inaccordance with the fire probability that takes into accountthe differences between the ease of fire occurrence andspread due to applications such as dining options,structure and residential buildings.
Since there is little knowledge about the evacuation in
the urban fire, elucidation of the complete behavior isdifficult. Therefore, in order to assess the risk with taking
into account the variety of actions the evacuees can take,the evacuation simulation model was constructed using
the approach of the agent-based modeling. Evacuationroutes of evacuees is basically the shortest path on theroad network towards the destination. However, evacuees
are also able to change their moving path in order to avoid
danger, depending on their perception of the.
We have carried out the evacuation simulation withdifferent fire-spread patterns, with several thousand trials
for each case.
3. CONCLUSIONS
Possibility for mass death caused by certain conditions
has been suggested.
Continuing from this, thesimulation results further reveal that certain places havehigher fatality rates, showing certain characteristics makea place a more dangerous for fire than other places.
Testing with different evacuation decision-making
models, consistency is visible among the results attained.
In the future, we hope to construct a more appropriate
evacuation model using experimental and questionnaires
and the like. We also plan to consider the sensitivity
analysis of the model.
REFERENCES[1] T. Kato, J. Yoshinaga, T. Shimura, T. Eda, Multi-times
simulation analysis of high risk during evacuation in a
post-earthquake urban fire (part1), presentation abstractsof Japan Association for Fire Science andEngineering,2015,pp.198-199.
[2]T. Shimura, T. Kato, J. Yoshinaga, T. Eda, Multi-times
simulation analysis of high risk during evacuation in a
post-earthquake urban fire (part2), presentation abstractsof Japan Association for Fire Science andEngineering,2015,pp.200-201.

― 45 ―

ANALYZING EFFECTIVE ROAD BLOCKING FOR IMPROVING EVACUATION TIMEFROM TSUNAMI USING TRFFIC FLOW MODEL
Ryosuke Yabe1, Yudai Honma2, Shigeki Toriumi3
1

Master Course Student, Department of Industrial and Management System Engineering, Waseda University, Tokyo,
Japan,
2
Lecturer, International Center for Urban Safety Engineering, The University of Tokyo, Tokyo, Japan,
3
Associate Professor, Information and System Engineering, Chuo University, Tokyo, Japan,
Correspond to Ryosuke Yabe (ryabe@toki.waseda.jp)

Keywords: Tsunami, Evacuation time, Road Closure, Traffic Jam, Bureau of Public Roads function
1. Introduction
In the Great East Japan Earthquake Disaster of March
11, 2011, a tsunami of unimaginable proportions assailed
the region, costing many lives. In its aftermath, tsunami
evacuation procedures came under review.
While there have been many studies of evacuation fromtsunamis, few studies have applied mathematical analysis
and few have analyzed traffic congestion during anevacuation.
This study attempts to reduce preferential car use in
major earthquake disasters by setting road closures. The
goal is to reduce the number of car users and detours
during emergency periods, leading to a reduction in traffic
congestion and evacuation time. The study is aimed atminimizing the total time taken for evacuation byproposing appropriate road closure strategies for a rangeof evacuation options.
2. EXISTING RESEARCH
Takada et al.[1] identified the impossibility of smoothevacuation using automobiles as one cause of increased
human casualties during a tsunami. They discussed thenecessity of post-analysis of evacuations by car whenplanning future tsunami evacuations strategies. Using theurban areas of Kesennuma City as the subject area, theyreconstructed the situation of the tsunami evacuation by
car during the Great East Japan Earthquake Disaster.
3. DESCRIPTION OF THE STUDY
We take Shizuoka City as our example, where tsunami
damage from a Nankai megathrust earthquake isanticipated. We used the basic road data available in the2005 edition of the National Digital Road Map Databasecreated by Sumitomo Electric as the network in this study.
Figure 1. Road network in the subject area

When a major earthquake strikes, evacuees will start toevacuate from a point in the red-lined area (start ofevacuation node) using roads (links) to evacuate to a pointin the black-lined area (high ground node). At the sametime, some of the roads will be closed to automobile trafficand evacuees are expected to recognize road closures inadvance.
If there are no road closures, a large percentage of
people will choose car, as that will allow them to evacuate
more quickly while fewer people will choose foot. For thisreason, road closures were introduced. These change theroutes available for car evacuation, producing a(perceived) increase in the time required. Under the roadclosure scenario, some evacuees who would havepreviously used a car will evacuate on foot or will take adetour, avoiding the concentration of traffic on certain
roads. This avoids major congestion, reducing the averagetime required to evacuate.
4. CALCULATION RESULTS
In this study, the average time required for evacuationin Shizuoka City improved from 9.41 min prior toapplying road closures to 6.62 min after application.
5.
FUTURE PLANS AND ANTICIPATEDCONCLUSIONS
These trial calculation results still leave some issues
unresolved. One is that the results are largely determinedby the parameters chosen. In particular, the number ofpeople and traffic capacity settings have a significant
impact and have not been appropriately allocated. In thisstudy, the number of vehicles in each link was calculated
as the cumulative total when deriving the number of carspassing through a link and calculating the average timerequired for evacuation. In reality, the flow of car traffic is
not even; hence, using a cumulative total will exaggerate
traffic congestion. A more accurate per minute calculation
will yield more accurate timings.
REFERENCES
Reference to a journal publication:[1] Takada, K., Fujiu, M., and Onomura, K., Research on
reproduction method of the situation of tsunami
evacuation with car – case study of Kesennuma city-.
Journal of Japan Society of Civil Engineers, A1(Structural and Seismic Engineering), Vol. 69, No.4, pp.
973-979, 2014.

― 46 ―

ÓßÌØÛÓßÌ×ÝßÔ ÓÑÜÛÔ ÚÑÎ ÛÏË×Ô×ÞÎ×ËÓ ßÒÜ ÑÐÌ×ÓËÔ ÔÑÝßÌ×ÑÒ ÑÚØÑËÍ×ÒÙÍ ßÒÜ ÖÑÞÍ É×ÌØ ÝÑÒÍÌÎß×ÒÛÜ ÝßÐßÝ×ÌÇÇò ØÑÒÓßï ô Çò ÓËÒÛÓßÍßî ô ¿²¼ Õò ×Óß×íïîÔ»½¬«®»®ô ×²­¬·¬«¬» ±º ×²¼«­¬®·¿´ Í½·»²½»ô Ì¸» Ë²·ª»®­·¬§ ±º Ì±µ§±ô Ì±µ§±ô Ö¿°¿²ß­­·­¬¿²¬ Ð®±º»­­±®ô Ü»°¬ò ±º ß®½¸ò Ý·ª·´ ú Û²ª·®±²³»²¬¿´ Û²¹òô Ì±µ§± Ü»²µ· Ë²·ª»®­·¬§ô Í¿·¬¿³¿ô Ö¿°¿²íÐ®±º»­­±®ô ×²­¬·¬«¬» ±º ×²¼«­¬®·¿´ Í½·»²½»ô Ì¸» Ë²·ª»®­·¬§ ±º Ì±µ§±ô Ì±µ§±ô Ö¿°¿²Ý±®®»­°±²¼ ¬± Ô»½¬«®»® Çò ØÑÒÓß ø§«¼¿·à··­ò«ó¬±µ§±ò¿½ò¶°÷Õ»§©±®¼­æ ¼·­¬®·¾«¬·±² ±º ¸±«­·²¹­ ¿²¼ ¶±¾­ô ¬®¿ª»´ ½±­¬ô ¼·­½±³º±®¬ ½±­¬
ïò×ÒÌÎÑÜËÝÌ×ÑÒ
×² ¬¸·­ °¿°»®ô ©» °®±°±­» ¿ ²»© ½¿´½«´¿¬·±² ³»¬¸±¼ ±º»¯«·´·¾®·«³ ¿²¼ ±°¬·³«³ ´±½¿¬·±² ±º ¸±«­·²¹­ ¿²¼ ¶±¾­¬¸»±®»¬·½¿´´§ ·² ¬¸» «®¾¿² ³±¼»´­ ©·¬¸ ½±²­¬®¿·²»¼ ½¿ó
°¿½·¬·»­ò É» ­¸±© ¬¸¿¬ ¿ °®±¾´»³ ©¸·½¸ ·²½±®°±®¿¬»­½±³³«¬·²¹ ¬®¿ª»´ ½±­¬ô ©±®µ·²¹ ¬®¿ª»´ ½±­¬ ¿²¼ ¼·­½±³óº±®¬ ½±­¬ ½¿«­»¼ ¾§ ¬¸» ¼»²­·¬§ ±º ¬¸» ²±¼» ©·´´ ¾» ¿²±²ó½±²ª»¨ ¯«¿¼®¿¬·½ °®±¹®¿³³·²¹ °®±¾´»³ò Þ§ «­·²¹
Î»º±®³«´¿¬·±²óÔ·²»¿®·¦¿¬·±² Ì»½¸²·¯«»ô ©» ¿®» ¿¾´» ¬±¼»®·ª» ¿²¿´§¬·½¿´´§ ¬¸» »¯«·´·¾®·«³ ¿²¼ ±°¬·³«³ ´±½¿¬·±²±º ¸±«­·²¹­ ¿²¼ ¶±¾­òËÎÞßÒ ÓÑÜÛÔÔ»¬ «­ ½±²­·¼»® ¿² «®¾¿² ½·¬§ ©¸·½¸ ½±²¬¿·²­ ±º ©±®µ·²¹°´¿½» øÉ÷ ½±²­·­¬·²¹ ±º ¬¸» ¾«·´¼·²¹ ²±¼»­ × øï · ×÷¿²¼ ®»­·¼»²½» °´¿½» øØ÷ ½±²­·­¬·²¹ ±º ¬¸» ¾«·´¼·²¹ ²±¼»­¶
Ö÷ô ¿²¼ ¿­­«³» ¬¸¿¬ ±²´§ ¬©± ¾«·´¼·²¹ ¬§°»­Ö øï»¨·­¬ ·² ¬¸» «®¾¿² ½·¬§ ·² ±®¼»® ¬± ³¿µ» º±®³«´¿¬·±² ³«½¸³±®» ½±²½·­»ò Ó±®»±ª»®ô °®±°±­¿´ «®¾¿² ³±¼»´ ·­ ·²¼»ó
°»²¼»²¬å ¬¸» ¬±¬¿´ ²«³¾»® ±º °»±°´» ³±ª·²¹ ÉóØ ·­ Ì øã½±²­¬¿²¬ ²«³¾»®÷ ¿²¼ ¬¸»®» ·­ ²± ·²A±© ¿²¼ ±«¬A±© ±º°»±°´» ¾»¬©»»² ±¬¸»® ½·¬·»­òÍ¿³» ¿­ ¬¸» °®»ª·±«­ ³±¼»´ ÅïÃô ©» ¿­­«³» ¬¸¿¬ ¬¸»®» ·­¬©± ¬®¿ª»´ ¬§°» ·² «®¾¿² ½·¬§ô ø·÷½±³³«¬·²¹ ¬®¿ª»´ ±º ÉóØ¿²¼ ø··÷©±®µ·²¹ ¬®¿ª»´ ±º ÉóÉŽ ¿²¼ ø···÷¼·­½±³º±®¬ ¼«» ¬±¼»²­·¬§ ±º ¬¸» ²±¼»ò ß¬ ¬¸·­ ¬·³»ô ±«® °«®°±­» ±º ¬¸·­ ­¬«¼§·­ ¬± ¼»®·ª» »¯«·´·¾®·«³ ¿²¼ ±°¬·³«³ ´±½¿¬·±² ·² ¬¸» «®¾¿²½·¬§ò Ñ² ¬¸¿¬ ¿½½±«²¬ôíò ÛÏË×Ô×ÞÎ×ËÓ ÔÑÝßÌ×ÑÒ
Ú·®­¬ô ©» ¼»®·ª» ¬¸» »¯«·´·¾®·«³ ¼·­¬®·¾«¬·±² ±º ¸±«­·²¹­¿²¼ ¶±¾­ò Ì¸·­ °®±¾´»³ ©·´´ ¾» ¿ ²±²½±²ª»¨ ¯«¿¼®¿¬·½ °®±ó¹®¿³³·²¹ °®±¾´»³ ­¬¿¬»¼ ¾»´´±©æ
³·²ò

Ö

Ì ã

¬¸·¶

øï

·

øï

¶

Ö÷ô

×ô ï

¶

Ö÷ô

øï÷

·ãï ¶ãï

ÚËÛ ã
×

×

×

î¼¸·¶ ¬¸·¶ õ·ãï ¶ãï

Ì ·ãï

· ãï

Ö

¸¶

©·

¦·© ø¨÷¼¨ õ

õ
·ãï ð

ï ©
¼ ©· ©·î ··¦¶¸ ø¨÷¼¨ øê÷
¶ãï ð

Ö

­ò¬ò

îò

×

Ö

×

¬¸·¶

øï

·

×÷

øé÷

¬¸·¶

øï

¶

Ö÷

øè÷

©· ã
¶ãï
×

¸¶ ã
·ãï

Ö

×

¬¸·¶

Ì ã

øç÷

·ãï ¶ãï

ð

©·

É·

øï

·

×÷

øïð÷

ð

¸¶

Ø¶

øï

¶

Ö÷

øïï÷

ð

¬¸·¶

øï

·

×ô ï

¶

Ö÷

øïî÷

Þ§ «­·²¹ Î»º±®³«´¿¬·±²óÔ·²»¿®·¦¿¬·±² Ì»½¸²·¯«» ÅîÃô ©»¿®» ¿¾´» ¬± ¼»®·ª» ¿²¿´§¬·½¿´´§ ¬¸» »¯«·´·¾®·«³ ´±½¿¬·±² ±º¸±«­·²¹­ ¿²¼ ¶±¾­òìò ÑÐÌ×ÓËÓ ÔÑÝßÌ×ÑÒÒ»¨¬ô ©» ¼»®·ª» ¬¸» ±°¬·³¿´ ¼·­¬®·¾«¬·±² ±º ¸±«­·²¹­ ¿²¼¶±¾­ò Ì± ­±´ª» ¬¸·­ °®±¾´»³ô ©» ±²´§ ¸¿ª» ¬± ½¸¿²¹» ¬¸»±¾¶»½¬·ª» º«²½¬·±² ¿­ º±´´±©­æ

Ö

¬¸·¶

©· ã

×

øî÷

¶ãï

×

¬¸·¶

øï

·

×÷ô

×

øí÷

õ

·ãï

·ãï

©¸»®» ¬¸·¶ ·­ ¿ ½±³³«¬·²¹ ²«³¾»® º®±³ · ¬± ¶ô ©· ·­ ¿©±®µ·²¹ ²«³¾»® ±º °»±°´» ¿¬ ²±¼» · ¿²¼ ¸¶ ·­ ¿ ®»­·¼»²½»²«³¾»® ±º °»±°´» ¿¬ ²±¼» ¶ò Ú«®¬¸»®³±®»ô ©» ¿­­«³» ¬¸»º±´´±©·²¹ ·²»¯«¿´·¬·»­ ¬± °®±ª·¼» ¬¸» ·¼»¿ ±º ½±²­¬®¿·²»¼½¿°¿½·¬§ô
ð
ð

©·
¸¶

É· ô
Ø¶ ô

©¸»®» É· ¿²¼ Ø¶ ·­ »¿½¸ ½¿°¿½·¬§ ¿¬ ¬¸» ²±¼»ò

øì÷
øë÷

×

¬¸·¶ ¼¸·¶ õ·ãï ¶ãï

×

¸¶ ã

Ö

ÚÍÑ ã î I©· ©··ãï · ãï

Ì

¼©
··

Ö

©·
¸¶
©· õ
¸¶
É·
Ø¶
¶ãï

øïí÷

ÎÛÚÛÎÛÒÝÛÍÅïÃ Ø±²³¿ô Çò ¿²¼ Õ«®·¬¿ô Ñòô îððìò Ü»¬»®³·²·²¹ Ö±·²¬Ü·­¬®·¾«¬·±²­ ±º Ø±³»­ ¿²¼ É±®µ°´¿½»­ Þ¿­»¼ ±² ¬¸»Î¿²¼±³ Ë¬·´·¬§ Ì¸»±®§ô Ö±«®²¿´ ±º ¬¸» Ý·¬§ Ð´¿²²·²¹×²­¬·¬«¬» ±º Ö¿°¿²ô íçô °°òéêçPééìò ø×² Ö¿°¿²»­»÷ÅîÃ Í¸»®¿´·ô Øò Üò ¿²¼ ß¼¿³­ô ÉòÐòô îðïíò ß Î»º±®³«´¿¬·±²ó
Ô·²»¿®·¦¿¬·±² Ì»½¸²·¯«» º±® Í±´ª·²¹ Ü·­½®»¬» ¿²¼ Ý±²ó¬·²«±«­ Ò±²½±²ª»¨ Ð®±¾´»³­ô Í°®·²¹»® ËÍò― 47 ―EVACUATION SIMULATION CONSIDERING BUILDING DAMAGES
S. KUBO1, R. FUJITA2, H. YOSHIDA3, M. L. L. WIJERATHNE 4, and M. HORI 5
Under Graduate Student, Faculty of Engineering, Kagawa University, Kagawa, Japan,
2
Maeda Corporation, Gifu, Japan,
3
Professor, Faculty of Engineering, Kagawa University, Kagawa, Japan,
4
Associate Professor, Earthquake Research Institute, The University of Tokyo, Tokyo, Japan,
5
Professor, Earthquake Research Institute, The University of Tokyo, Tokyo, Japan,
Keywords: Integrated Earthquake Simulator, Evacuation simulation, Structure response analysis, Building damages
1

1. INTRODUCTION
Serious damages were occurred by the Great East Japan
Earthquake in March 11th, 2011. Most victims were
engulfed by Tsunami, which were caused by the blindness
in the decision of evacuation and the road block bycollapsed buildings and stalled vehicles. In order to solve
the evacuation problems, ‘Integrated Earthquake
Simulator’ (IES) is developed by the Earthquake ResearchInstitute of the University of Tokyo. The IES enables us todiscuss the circumstances or state of city after anearthquake comprehensively through the seamless
analyses. In this study, the evacuation analysis isconducted by the IES to discuss the dependency ofevacuation on the existence or non-existence of road block
collapsed buildings.
2. STUDY FLOWIn this study, first, a structure response to observed
earthquake waves (the acceleration waveform at groundlevel) for the city model of Takamatsu is analyzed by the
IES. The numerical result is the time history displacementof every floor of each structure. From the result story drift(ߠ) is calculated by the following equation.
ߠ ൌ ߜȀ‫ܪ‬where ߜ and H are the horizontal displacement andheight of each layer, respectively.
Secondary, the city model in which building damages
are considered is prepared and evacuation simulations
considering building damages caused by the earthquake
are carried out. The numerical results show that theevacuation time depends on the existence or non-existenceof damages and the difference extent of damages. Thus,
the area in which there is no route for evacuation is set up
with taking the damage area into account (see Figure 1 andFigure 2). The damage area is determined by the story drift.
For example, a damage is assumed to be occuered if thestory drift exceed a certain value.

Figure 1. City model ignoring building damage

Figure 2. City model concerning building damages
3. ANALYSIS RESULTS
In the-area where low, middle and high-rise buildings
are mixed, the building damages concentrate on the
middle rise buildings when the Southern Hyogo prefectureearthquake is employed as an incident seismic wave. From
the result, it is clarified that the road at the periphery ofthe damaged buildings can be blocked or narrowed.
Meanwhile, the building damages concentrate on the
low-rise buildings when the Great East Japan Earthquake
is employed. Since the road around the low-rise buildings
is generally narrow, some of the agents go by aroundabout evacuation route or change the evacuation
place.
In the area where the low-rise buildings is dominant, thebuilding damages are limited to a portion of middle-risebuildings when the Southern Hyogo prefecture earthquake
is employed. Therefore, the road is rarely blocked or
narrowed by the building damages.
On the other hand, the road block or narrowing by thebuilding damages is confirmed everywhere when the
Great East Japan Earthquake is employed. Due to the roadblock or narrowing, the agents go by a roundabout
evacuation route and head back. At time like this, it issuggested that it may take much time to evacuate.
4. CONCLUSIONS
The analysis results suggest that the evacuation action
can depend on the damage situation of buildings, and thatthe introduction of the visualization of damaged buildings
enables us easily to grasp the damage situation. Thus, theproposed scheme is effective to discuss and evaluate thenumerical results.
REFERENCES[1] M. Hori, H. Miyajima, Y.Inukai, K. Oguni, Agent
Simulation for Prediction of Post-Earthquake Mass
Evacuation, Japan Society of Civil Engineers, A, Vol.64,
No.4, pp.1017-1036, 2008.

― 48 ―

A 2D-3D HYBRID MODEL BASED ON STABILIZED FINITE ELEMENT METHOD FORTSUNAMI RUNUP SIMULATION
G. LING1, K. KASHIYAMA2, and J. MATSUMOTO3
Ph.D. Student, Department of Civil and Environmental Engineering, Chuo University, Tokyo, Japan,
2
Professor, Department of Civil and Environmental Engineering, Chuo University, Tokyo, Japan,
3
Senior Researcher, National Institute of Advanced Industrial Science and Technology (AIST), Ibaraki, Japan
Correspond to Prof. K. KASHIYAMA (kaz@civil.chuo-u.ac.jp)
1

Keywords: 2D-3D Hybrid Model, Stabilized Finite Element Method, VOF Method, Tsunami
1. INTRODUCTION
In recent years, a number of destructive tsunamis reveal
the importance of making timely predictions of tsunami
runup and giving tsunami warnings. Since the traditional
2D tsunami simulations which are based on thehydrostatic approximation may be limited to get the detail
information around structures on the coastal area and the3D tsunami simulations may cost too much, we propose a
2D-3D hybrid model to compute the wave propagation by
2D and the area around structures by 3D. The present
method is applied to several benchmark examples to show
the validity and efficiency of the method.

wave propagation at 0.6s by using the 2D-3D hybrid
model. Figure 4 shows a comparison of the 2D-3D hybrid
models, 2D models and the experimental results.
According to Figure 4, we can see the results of the 2D3D hybrid model show better agreement with theexperimental results.

2. NUMERICAL METHODSThe 2D-3D hybrid model couples the 2D shallowwater/Boussinesq equation with the 3D Navier-Stokes
equation are based on the stabilized finite element method[1]. The advection equation is applied to the interfacecapturing approach that based on the VOF method [2]with the 3D analysis. The linear triangular element is
employed for the spacial discretization, the CrankNicolson method is employed for the temporal
discretization, and the Element-by-Element Bi-CGSTAB
method is employed for solving the simultaneous linear
equations.
A 2D-3D overlapping method which has considered
about the conservation of the flow rate is employed.
Figure 1 shows the images about the coupling domain.
The velocity and the water depth which are computed by
the 2D domain will be used as the 3D inflow boundary
conditions. As the same, the velocity and the water depth
which are computed by the 3D domain will be used as the2D outflow boundary conditions.

Figure 2. The runup of solitary waves problem

Figure 1. 2D-3D overlapping method
3. NUMEIRCAL EXAMPLES
In order to evaluate the effectiveness of the 2D-3D
hybrid model, we applied the runup of solitary waves
problem [3] shown in Figure 2. We set the solitary wave
with initial amplitude
climbing over thebeach with. Figure 3 is a snapshot of the

Figure 3. A snapshot of the wave propagation (0.6s)

Figure 4. Surface elevation profiles
3. CONCLUSIONSIn this paper, a 2D-3D hybrid model based on thestabilized finite element method for tsunami runup
simulation is presented. By comparing the numericalresults with the experimental results, it can be concluded
that the present method is useful for the tsunami runup
simulation. Our future work is the application to real city
terrains.
REFERENCES[1] S. Takase, K. Kashiyama, S. Tanaka and T. Tezduyar,
Space-time SUPG formulation of the shallow waterequations, Int. J. Numer. Meth. Fluids 64 (2010) 1379–
1394.
[2] C.W. Hirt, B.D. Nichols, Volume of fluid method forthe dynamics of free boundaries, J. Comput. Phys. 39(1981) 201-225.
[3] C.E. Synolakis, The runup of solitary wave, J. Fluid
Mech. 185 (1987) 523-545.

― 49 ―

APPLICATION OF THE FEMA P-58 METHOD FOR REGIONAL SEISMIC LOSS PREDICTION
X.Z. LU1, X. ZENG2, T.Y. YANG3, and Z. XU4
1

Professor, Department of Civil Engineering, Tsinghua University, Beijing, P.R. China,
Ph.D student, Department of Civil Engineering, Tsinghua University, Beijing, P.R. China,
3
Associate Professor, Department of Civil Engineering, University of British Columbia, Vancouver, Canada,
4
Associate Professor, Sch. Civil & Environm. Engn., University of Science and Technology Beijing, Beijing, P.R. China,
Correspond to Prof. Xinzheng LU (luxz@ tsinghua.edu.cn)
2

Keywords: earthquake engineering; FEMA P-58; earthquake economic loss; regional seismic damage simulation

2. LOSS PREDICTION METHODOLOGYIn FEMA P-58 method, the building seismic loss is thesum of the repair cost of all components. The repair costof each component is calculated through the following
procedure: (a) Obtain the EDPs of the building story inwhich the component is located; (b) Determine thedamage state utilizing the fragility curves assigned to thecomponent; (c) Determine the repair cost according to
the total quantity and unit repair cost of the component.
The EDPs are obtained through nonlinear time history
analysis (THA) using the MCS model. The parameters of
MCS models can be conveniently determined for a large
number of buildings with different structural types[2].
Moreover, the THA of a region using the MCS model
can be very efficient[3]. Hence the adoption of the MCS
model overcomes the difficulty of EDP calculation.
3. CASE STUDY
A case study of Tsinghua University campus in Beijing
is performed to demonstrate the application of FEMA P58 method for regional seismic loss prediction. The
region consists of 619 buildings. Masonry structure
comprise 55%, others are reinforced concrete frame andshear wall structures. The seismic loss prediction is
performed at three earthquake hazard levels. Figure 1
shows the result. The total seismic loss ratio is defined as
the ratio of total seismic loss to the total replacement

value (US$ 7.476 billion) of the region. More details
such as the loss results of components can be predictedby the proposed method for each building (Figure 2).
1

Probability of non-exceedance

0.75
Service Level
Earthquake (SLE)

0.5

Design Base
Earthquake (DBE)
0.25

Maximum Considered
Earthquake (MCE)

0
0%

10%

20%

30%

40%

50%

Total sesmic loss ratio

Figure 1. Cumulative probability of the total earthquake
loss ratio
0.6%

Loss ratio

1. INTRODUCTIONFEMA P-58 is the state of the art of seismic loss
prediction method[1]. Compared to the existing buildingbased loss prediction method such as HAZUS, the
FEMA P-58 method directly take the fragility of every
structural and nonstructural component in a building into
account, hence it can generate a more realistic loss result.
The FEMA P-58 method has been successfully applied
to individual buildings. However, currently there is nosuch application for regional predictions. Therefore, apracticable approach for regional seismic loss prediction
based on the FEMA P-58 method is proposed in thisstudy. The key challenge when using the FEMA P-58
method for regional seismic loss predictions is thedifficulty of obtaining the engineering demandparameters, or EDPs (i.e. building seismic responses,
such as inter-story drift ratio or floor acceleration) for alarge number of buildings within the region. In this study,
the multi-story concentrated-mass shear (MCS) model[2]is adopted for efficient seismic response analysis.

4th story
3rd story

0.4%

2nd story
1st story

0.2%
0.0%

Figure 2. The loss details of nonstructural PGs at eachstory of an individual building at the MCE hazard level.
4. CONCLUSIONS
A component-based regional seismic loss prediction
method is proposed based on the FEMA P-58 method
and the MCS model. The outcomes of this work can beused as a reference for urban seismic disaster mitigationor earthquake insurance.
REFERENCES[1] Federal Emergency Management Agency, Seismic
Performance Assessment of Buildings, Tech. Rep. FEMA
P-58, Washington, D.C., 2012.
[2] X.Z. Lu, B. Han, M. Hori, C. Xiong, Z. Xu, A
Coarse-Grained Parallel Approach for Seismic Damage
Simulations of Urban Areas Based on Refined Models
and GPU/CPU Cooperative Computing. Adv. Eng. Softw.
70 (2014) 90-103.
[3] Z. Xu, X.Z. Lu, H. Guan, B. Han, A.Z. Ren, Seismic
Damage Simulation in Urban Areas Based on a HighFidelity Structural Model and a Physics Engine. Nat.
Hazards 71(2014), 1679-1693.

― 50 ―

CODE2015

Technical Program – Day 2

Session 4

DEVELOPMENT OF A VOLUNTARY DECISION-MAKING MODEL FOR RESIDENTS DURING
POST-DISASTER RECOVERY
Yasmin BHATTACHARYA1 and Takaaki KATO2
1
2

Project Researcher, Institute of Industrial Science, University of Tokyo, Tokyo, Japan,
Associate Professor, Institute of Industrial Science, University of Tokyo, Tokyo, Japan,

Keywords: Agent-based simulation, disaster recovery, decision-making model
1. INTRODUCTION
The process of disaster recovery involves the dynamic
interaction of multiple factors which often remains hiddenfrom the surface[1]. Furthermore, due to empiricallimitations, recovery outcomes are often superficially
understood.
One of the major sources of the dynamics in the system
is residents’ behavioral decisions. The way affected
residents make decisions under uncertainty and how thesedecisions get propagated over space and time toeventually affect the overall recovery in the community isnecessary to be understood in order to have an effective
recovery policy. Therefore, the objective of this researchis to introduce a new methodological approach to capture
the dynamics of post-disaster recovery by modelingresidents’ decision making and interactions within therecovery environment.
2. MODEL COMPONENTS
This paper presents additional enhancement on a multiagent simulation model for recovery[2], which isdeveloped by taking into consideration the pre-disastercommunity characteristics of a disaster-hit region; thetype of damage; and the recovery policies applied to thesituation. In this approach, affected households aremodeled as agents with certain requirements and priorities,
and a belief structure which is updated as the recovery
proceeds. The interaction between their requirements andbelief forms the basis for their decision to move awayfrom or stay in the recovery area.
By using a multi-agent model, macro level
phenomenon can be understood in terms of micro level
decision-making of agents. In this research we haveparticularly focused on the decision-making inconsideration to the voluntary permanent relocation of thehousehold from the temporary housing and not considered
the processes involved in the immediate relocation ofhouseholds from the disaster area to temporary locations.
The decision-making for permanent relocation involves
different kinds of characteristics of the agent which can bebroadly divided into Primary and Secondary agent
characteristics. Primary characteristics relate to the levelof attachment to a place; minimum infrastructure
threshold; and affordable cost range of an agent. While theSecondary characteristics relates to the preferences of thehousehold regarding the different attributes of an areasuch as community, employment, and the existence ofschools and hospitals etc. which act as weightages whencalculating a household’s perception of attractiveness fora particular area.

Using this value of attractiveness as a measure, theagents can decide whether to remain in the temporary
housing and wait, move to the external area, or return tothe disaster area at every one month interval looking at the
progress of recovery.
A dynamic component in the model is the decline of the
Ideal Attractiveness level which is the total area
attractiveness value the agent hopes its final location will
have. This ideal begins to declines as the agent continues
to be frustrated in the temporary housing unable to find apermanent solution.
Finally, Recovery Policy measures are also represented
in the model in terms of short-term and long-term,
corresponding to the immediate response measuresprovided after the disaster and the permanent recovery
measures provided, respectively. Short-term measures
include temporary housing locations, duration of stay andquality of the area in term of area attributes discussed
previously. Long-term measures include the sharing/nonsharing of the reconstruction plan (including the final
quality of recovery and the intended date of completion)with the citizens, employment policies and subsidies forrebuilding given to the victims.
3. CONCLUSIONS
Through running the simulation, many of theempirically observed phenomena were validated.
Furthermore, it was made apparent that even if aphenomenon seems random at one level of focus, it does
not mean that there is no influencing factor on a different
level of focus. Since this simulation enables the attainmentof a micro-level data in time-sequence it can further reveal
the difference in movement patterns between different
types of population depending on a recovery policy evenwhen the final outcome is same. Additional simulations
may yield insights about the nature of recovery that are
not empirically observable.
REFERENCES[1] A. Nejat and I. Damnjanovic, Agent-Based Modelingof Behavioral Housing Recovery Following Disasters,
Computer-Aided Civil and Infrastructure Engineering 27,
no. 10 (2012) 748–763.
[2] Y. Bhattacharya and T. Kato, Deciphering Recovery
Characteristics through a Multi-agent Simulation Model,
Proceedings of 3rd International Conference on Urban
Disaster Reduction, Earthquake Engineering ResearchInstitute, Boulder, CO, USA, 2014.

― 51 ―

STABILIZED FINITE ELEMENT METHOD BASED ON VOF METHODFOR FREE SURFACE FLOW USING LES
M. OTA1, G. LING2, and K. KASHIYAMA3
1

Graduate Student, Department of Civil and Environmental Engineering, Chuo University, Tokyo, Japan,
2
Ph.D. Student, Department of Civil and Environmental Engineering, Chuo University, Tokyo, Japan,
3
Professor, Department of Civil and Environmental Engineering, Chuo University, Tokyo, Japan,
Correspond to Prof. K. KASHIYAMA (kaz@civil.chuo-u.ac.jp)

Keywords: Stabilized Finite Element Method, VOF Method, LES, Fluid Force, Navier-Stokes Equation
1. INTRODUCTION
Tsunami causes the serious damage to the structures incoastal area. In order to compute the complicated behaviorof tsunami around coastal structure accurately, it isnecessary to employ the three-dimensional equations
based on Navier-Stokes equation and turbulence model.
The finite element method is a useful discretization
method for the arbitrary computational domain.
This paper presents a stabilized finite element method
for the free surface flow using the LES based on the
Smagorinsky model [1]. The present method is applied tothe three-dimensional dam break problem with the squarepillar.
2. NUMERICAL METHOD
The filtered incompressible Navier-Stokes equation andcontinuity equation are employed for the governing
equations of flow field. The stabilized finite element
method based on the SUPG/PSPG method [2] is employed
for the discretization in space. The implicit method basedon the Crank-Nicolson method is employed for thediscretization method in time. The fluid force acting on
the structure is evaluated by the weak form of thegoverning equations.
The VOF method [3], which is one of the interfacecapturing methods, is employed as the free surface flow
model. The interface function is governed by the timedependent advection equation. The SUPG method is
employed for the discretization in space, and the CrankNicolson method is employed for the discretization
method in time.
3. NUMERICAL EXAMPLE
The present method is applied to the three-dimensionaldam break problem with the square pillar [4]. The
computational model and initial condition are shown in
Figure 1 (left). The mesh with 7,029,070 tetrahedron
elements and 1,224,746 nodes is used in this problem. The
minimum mesh size and time increment are 1.27×10-3 m
and 0.001 s. The Smagorinsky constant is assumed to be
0.10.
Figure 1 (right) shows the profile of free surface and thevalues of pressure on the structure at t=0.40 s in the caseusing LES and slip boundary condition. The color denotes
the distribution of pressure.
Figure 2 shows the time history of the drag force acting
on the structure. From this figure, the result using the LES

is good agreement with the experimental results. Also, thecomputed results using the slip boundary condition gives
a better solution comparing with the non-slip boundary
condition.

Figure 1. The computational model (left) and the profile
of free surface and the values of pressure on the structure
at t=0.40 s (right)

Figure 2. Comparison of the time histories of drag force
acting on the structure
4. CONCLUSIONSThe numerical method for the free surface flow using
LES based on the stabilized finite element method has
been presented. From the results, the present method isuseful for the free surface flow analysis.
The development of the fluid-structure interaction
method is left for the future work.
REFERENCES[1] J. Smagorinsky, Monthly Weather Review, 91, 3, 1963,
pp. 99-164.
[2] T.E. Tezduyar, Advances in Applied Mechanics, 28,
1992, pp. 1-44.
[3] C.W. Hirt, B.D. Nichols, Journal of Computational
Physics, 39, 1981, pp. 201-225.
[4] M. Gomez-Gesteira, R.A. Dalrymple, Journal of
Waterway, Port, Coastal and Ocean Engineering, 130,
2004, pp. 63-69.

― 52 ―

PRE-ESTIMATION METHOD FOR LANDSLIDE HAZARD WITH MONITORING ANDSIMULATION
J.P. ZHANG1, 2, J.G. CHEN1, 2 and H.Y. YUAN1, 2
Institute of Public Safety Research, Tsinghua University, Beijing, China,
Department of Engineering Physics, Tsinghua University, Beijing, China,
Correspond to Dr. J.P. ZHANG (jp_zhang@tsinghua.edu.cn)

1
2

Keywords: Landslide hazard, Pre-estimation, Monitoring, Simulation
1. INTRODUCTION OF THE PRE-ESTIMATIONMETHOD
Although it is very difficult to predict the exact time ofa landslide hazard, it doesn't mean we can do nothing
about it. If one can obtain enough information on what
will happen or the impact degree in possible landslide
hazards, corresponding protective measures can be taken
previously to avoid or reduce losses.
The pre-estimation process can be realized nowadays
with the help of the modern information technology. The
basic techniques include in-situ monitoring[1] andcomputer simulation. Monitoring is an important way to
get real time data from the slope concerned, with sensors
and network transmission. The monitoring data include
deformation, displacement, cracks, vibration and theinducing variables of the environment, e.g. rain fall, floodor earthquake. These data are the foundation of the preestimation process, which are quite essential for theaccumulation of the database, the verification andmodification of the analytical model and the evaluation
for the current state of the slope. Computer simulation isbased on the big data collected through monitoring andhigh performance computing. Usually, the analysis
method is based on empirical models established byanalyzing the data collected to build the relationshipbetween the state of the slope and the monitoring data.
However, more accurate and effective pre-estimation
should combine with numerical simulation[2]. With thedevelopment of computer technology and modern
computational mechanics[3], the deformation, fracture andmovement of the slope under various external loads canbe simulated vividly with reasonable and objective results
if the right inputs and parameters are given. The
monitoring data can be used to validate the computational
method and provide information to determine the current
state and necessary parameters of the slope for thesimulation. Thus, the failure process and damage scale ofthe landslide hazard under distinct possible conditions canbe obtained directly through monitoring data andcomputer simulation.
2. IMPLEMENTATION SCHEME OF THE PREESTIMATION METHOD
The combination of in-situ monitoring and numericalsimulation is the key feature of the pre-estimation method.
The implementation scheme is shown in figure 1. The
monitoring data provides basic information(e.g., initialand boundary conditions) for simulation and verifies thenumerical model. Meanwhile, the numerical model

determines the most essential parameters to monitor andoptimizes the monitoring scheme specifically. Then, one
can preview possible landslide hazards under distinct
given conditions and pre-estimate the possible damage
degree. Mutual verification and modification ofmonitoring scheme and numerical model ensures therationality of this method.

Figure 1. Scheme of the pre-estimation method
3. CONCLUSIONSThe pre-estimation method for landslide hazards within-situ monitoring and numerical simulation can make fulluse of the monitoring data and the latest computing
technology. Preview of possible landslide hazards in given
conditions can be called "condition prediction". It can
provide suggestions for the design of the anti-landslide
structures and help government to make proper plans anddecisions in landslide-prone regions.
REFERENCES[1] X. J. Yang, et al. Stability and remote real-timemonitoring of the slope slide body in the Luoshan mining
area. International Journal of Mining Science andTechnology. 25(2015)761-765.
[2] Z. Li, et al. A case study integrating numericalsimulation and GB-InSAR monitoring to analyze flexural
toppling of an anti-dip slope in Fushun open pit.
Engineering Geology, 197( 2015) 20-32.
[3] T. Łodygowski, J. Rakowski, P.Litewka. Recent
Advances in Computational Mechanics. CRC Press,
Taylor & Francis Group, 2014.

― 53 ―

PREPARING SAFE URBAN ROADS FOR FLOOD AND EARTHQUAKE IN YANGON,
MYANMAR
T. SRITARAPIPAT1, W. TAKEUCHI2
Ph.D. Student, Institute of Industrial Science, The University of Tokyo, Tokyo, Japan,
Associate professor, Institute of Industrial Science, The University of Tokyo, Tokyo, Japan,
Correspond to T. SRITARAPIPAT (tanakorn@iis.u-tokyo.ac.jp)
1

2

Keywords: Road, Flood, Earthquake, Traffic Congestion
1. INTRODUCTION AND METHODOLOGY
Around the word, disaster is the serious problem thatnegatively impacts to humanity. Especially, it happens inan urban area that has the high density of population anda number of human activities. Road network plays animportant role in transportation for people, goods,
objects and etc.. Also, the urban road system is rather
crucial as a media to link all the human activities. Whenroads are rapidly damaged by the disaster and not able tocommonly use, it has a serious effect on transportationand possibly breaks down some related activities. In
order to reduce effects of disasters such as flood andearthquake in the urban area, it is necessary to prepare
safe roads which are not damaged by the disaster tosecurely evacuate or to normally operate when disasteroccurring.
Estimating hazard roads for the flood was introduced
by Kalantari [1]. They used the factors from soil wetness
index, road density, soil properties and channel slope toevaluate the hazard roads. Preparing pedestrians toevacuate in case of the earthquake was proposed by
Bernardini [2]. They applied previous databases onhuman evacuation behaviours to find safe pedestrians.
Moreover, traffic state estimation was introduced by
Shan [3]. They employed local camera with big GPS
data to detect traffic congestion.
This work proposed a methodology to extract safe
roads for flood and earthquake in Yangon, Myanmar. In
this research, Stereo GeoEye Images and Landsat timeseries were used to provide flood and earthquake maps[4]. While roads were manually extracted by ICUS
(International Center for Urban Safety Engineering). The
roads with disaster risk were computed by integrating
between road products and disaster risk map.
For the factors to indicate flood risk, there are two
factors; (1) land elevation, (2) floodway. For the factors
to indicate earthquake risk, there are three factors; (1)land slope, (2) indirect old building, (3) distance fromfaultine. Moreover, traffic congestion can be indirectly
estimated by observing land value area. There are two
factors; (1) urban expansion (2) the height of building toindicate land value. After processing for disaster risk
map, we integrated the roads with a buffer zone anddisaster risk map to produce the map of the roads withflood and earthquake risks.
2. EXPERIMENTAL RESULT
In our experimental results, there are two safe roadmaps for flood and earthquake in figure 1 and 2. Also,
indirect traffic congestion map was illustrated in figure 3.

(a)

safe

(c) low traffic

danger(b) safe

dangerhigh traffic

Figure 1. Safe Roads for (a) flood, (b) Earthquake and (c)traffic state
According to the resultant maps, the roads, which arelocated in high elevations and non-floodway, are secure
for the flood. The roads which are located in highelevations and non-floodway are secure when earthquake
occurring. Additionally, the roads, which are near thecommercial building and are the center of the urban area,
are high traffic state. Combining dangerous roads fromdisaster with traffic state, the huge impact from disasterin transportation can be detected.
3. CONCLUSIONSThis research introduced a method to estimate safe
roads for flood and earthquake in Yangon, Myanmar.
Satellite images from GeoEye and Landsat were
employed to provide flood and earthquake maps.
Integrating between roads with a buffer zone and disaster
map, safe road maps for flood and earthquake were
proposed. These maps can be used to prepare andmitigate the effect of the disaster in transportation.
REFERENCES[1] Z. Kalantari, A. Nickman, S. W. Lyon, B. Olofsson, L.
Folkeson, J. of Env. Management, 133 (2014) 69–77.
[2] G. Bernardini, E. Quagliarini, M. D'Orazio, Safety
Science, 82 (2016) 77–94.
[3] Z. Shan, Q. Zhu, Neurocomputing, 169, 2 (2015)
134–143.
[4] T. Sritarapipat, W. Takeuchi, ACRS2015, October
2015.

― 54 ―PROBABILISTIC ANALYSIS OF TSUNAMI HAZARDUSING NUMERICAL SIMULATIONS AND RESPONSE SURFACE
T. KOTANI1, S. TAKASE2, S. MORIGUCHI3, K. TERADA4, Y. FUKUTANI5, Y. Otake6, K. NOJIMA7,
M. SAKURABA7
Graduated student (M2), Graduate school of Engineering, Tohoku University, Sendai, Japan,
2
Assistant Professor, School of Engineering, Tohoku University, Sendai, Japan,
3
Associate Professor, International Research Institute of Disaster Science, Tohoku University, Sendai, Japan,
4
Professor, International Research Institute of Disaster Science, Tohoku University, Sendai, Japan,
5
Tokio Marine & Nichido Risk Consulting Co., Ltd., Chiyoda, Japan,
6
Associate Professor, School of Engineering, Niigata University, Niigata, Japan,
7
Research and Development Center, Nippon Koei Co., Ltd., Tsukuba, Japan,
1

Keywords: Tsunami Hazard analysis, Numerical Simulation, Response Surface
heights were calculated. Consequently, we infer that theresult is consistent against observed data.
Uncertainty analysis

Empirical
knowledge on
uncertainty by database

Analytical modeland parameters
Tsunami analysis

2. OUTLINE OF THE PROPOSED METHOD
A framework of the reliability based design forgeostructures had been proposed by Honjo [2]. In thisstudy, the basic concept of the framework is applied totsunami hazard evaluation. Fig.1 shows a flow of the
proposed method. As shown in the figure, evaluation
procedures are separated into three parts; (i) numericalsimulation, (ii) uncertainty analysis, and (iii) reliability
assessment. In the part of the numerical simulation, aseries of simulations are performed under the different
calculation conditions to obtain the response surface(RS) which is a relation between evaluation target andbasic variables. Then, in the part of the uncertainty
analysis, the uncertainty of the basic variables arequantified based on the database and engineeringknowledge. The Monte Carlo simulation (MCS) is
finally performed using RS and quantified uncertainty ofthe basic variables to obtain the probability density
distribution of evaluation target.
3. APPLICATION EXAMPLE
The proposed method was applied to evaluate thetsunami that had been induced by the great east Japan
earthquake. Fig.2 shows obtained probability density
distribution of maximum tsunami height at Sendai and
Ishinomaki. Observed maximum tsunami heights at eachpoint obtained from the results of the field survey arealso shown in the figure. In addition, the exceedance
probability against the observed maximum tsunami

Basic variables ‫ݔ‬Tsunamisimulation- Tsunami analysis toolܻ ൌ ‫ܨ‬ሺ‫ݔ‬ሻ
- Engineering judgment

Quantification of
uncertainty in ‫ݔ‬Uncertainty
analysis

1. INTRODUCTION
A method of probabilistic tsunami hazard analysis had
been proposed by the Japan Society of Civil Engineering(JSCE) [1], and the method have been widely used insubsequent studies. In the method, however, numericalsimulations are only partly used to evaluate theprobability density function. By using numericalsimulations more efficiently, more accurate method canbe potentially established. The purpose of this study is to
propose a new approach which can enable us to use
numerical simulations properly in the probabilistic
tsunami hazard analysis.

ex.
- Tsunami height- Time of tsunami-arrival- To be flooded or not

Response of evaluation
phenomenon : ܻ

Reliability
assessment

Random number
generation: ‫ݔ‬

Reliability
assessment(MCS)

Unsertainty
assesmentof ܻ:
ex. ܲ ܻ ൐ ܻ ‫ כ‬ൌ ܲ௙

Calculating
Exceedance probability

Approximationby RSܻ ൎ ‫ܨ‬ሺ‫ݔ‬ሻ

Figure 1. A flow of proposed method

Field survey result(6.902m)

Field survey result(10.833m)

60.6%

54.6%

Figure 2. Probability density distribution of maximum
tsunami height at Sendai (left) and Ishinomaki (right)
4. CONCLUSION
The new method of probabilistic tsunami hazard
analysis was proposed, and the method was applied toevaluation of the tsunami in the great east Japan
earthquake. Although consistent result was obtained, it isnecessary to perform more detailed validation to clarify
the capability and the advantage of the proposed method.
REFERENCES[1] The Tsunami Evaluation Subcommittee, The Nuclear
Civil Engineering Committee, JSCE, “Method of probabilistic tsunami hazard analysis”, 2002. (in Japanese)[2] Honjo, Y., “Challenges in Geotechnical Reliability
Based Design”, proceedings of 3rd International Symposium on Geotechnical Safety and Risk, pp.11-27, 2011.

― 55 ―

TOWARDS AN EXPECTATION OF TSUNAMI FEATURES USING SAR DATA
H. GOKON1, S. HAYASHI2, S. KOSHIMURA3 and K. MEGURO4
1

Assistant Professor, Institute of Industrial Science, The University of Tokyo, Tokyo, Japan,
2
Graduate Student, Graduate School of Engineering, Tohoku University, Sendai, Japan,
3
Professor, International Research Institute of Disaster Science (IRIDeS), Tohoku University, Sendai, Japan,
4
Professor, Institute of Industrial Science, The University of Tokyo, Tokyo, Japan,
Correspond to Assistant Professor H. Gokon (gokon@iis.u-tokyo.ac.jp)
Keywords: Tsunami Feature, Remote Sensing, Synthetic Aperture Radar, Tsunami Fragility Function
1. Introduction
Field survey after an extensive tsunami disaster needsa large amount of effort and times. Furthermore, thedelay of the field survey might cause the lost of theevidence of tsunami features. To solve this problem, thisstudy aims at developing a technology to expect thetsunami features in a tsunami affected area using remotesensing technology.
1. Study area and data set
The study area encompassed Arahama, Sendai ciy in
Miyagi prefecture, which was devastated by the 2011
Tohoku earthquake and tsunami event. Pre- and postevent high-resolution Synthetic Aperture Radar(TerraSAR-X) data, that observed affected areas on 20
October 2010 (UTC) and 12 March 2011 (UTC), wereused. Built-up areas were identified by a GIS data withbuilding footprint published by Zenrin company. The
result of the expectation of tsunami features was verified
by the spatial distribution of tsunami in- undation withmaximum inundation depth and the maximum flow
velocity, that were calculated by Hayashi et al.(2013)[2].
2. METHOD
Expectation of tsunami features consists of 1) Building
damage estimation with TerraSAR-X data, 2)
Conversion of the estimated building damage to thefeatures of tsunami inundation using tsunami fragility
function.
2.1.BUILDING DAMAGE ESTIMATIONFirst, damage estimation of washed-away buildings is
conducted by applying the following 6 procedures to preand post-event TerraSAR-X data ([1]); 1) Pre-processing,2) Change detection, 3) Extraction of built-up areas, 4)
Object-based image processing, 5) Estimating damage
ratio of washed-away buildings.
2.2.Conversion of damage ratio to tsunami features
Based on the fragility function proposed by Hayashi etal.(2013)[3], the damage ratio estimated by the damage
function was changed into maximum inundation depth
and maximum flow velocity. The fragility functions forwooden buildings, that were used for estimating tsunami
features, are shown in. Finally, the expected values of theFig1. An example of expecting flow depth using SAR
data (A) Expected, (B) Calculated
flow depth and flow velocity was verified by the results
proposed by Hayashi et al.(2013)[2].
3. RESULTS AND CONCLUSIONSThe summary of this study is shown as follow; (1) A
new method was proposed to expect the maximum
inundation depth and maximum flow velocity, by using
pre- and post- event TerraSAR-X data, and the feasibility
was tested. (2) It could be found that, as the damage ratio
is close to one, the overestimation was easy to be caused.
(3) On the other area, the tsunami features could be
estimated almost well.
REFERENCES[1] H. Gokon, S. Koshimura and M. Matsuoka, An
Object-oriented Method to Detect Building Damage
using TerraSAR-X Data, Journal of Japan Society of
Civil Engineers, Ser. B2 (Coastal Engineering), Vol.70,
No.2, 2014, pp.1481-1485.
[2] S. Hayashi and S. Koshimura, The 2011 Tohoku
Tsunami Flow Velocity Estimation by the Aerial Video
Analysis and Numerical Modeling, Journal of Disaster
Research, Vol.8, No.4, 2013, pp. 561-572.
[3] S. Hayashi, Y. Narita and S. Koshimura, Developing
Tsunami Fragility Curves from the Surveyed Data andNumerical Modeling of the 2011 Tohoku Earthquake
Tsunami, Journal of Japan Society of Civil Engineers,
Ser. B2 (Coastal Engineering), Vol.69, No.2, 2013,
pp.386-390.

― 56 ―

CODE2015

Technical Program – Day 2

Session 5

NUMERICAL SIMULATION OF HIGH VELOCITY IMPACT ON STEEL-PLATE CONCRETE
WALL
J.W. JUNG1, H.W. JANG2, and J.W. HONG3
Ph.D. Candidate, Civil and Environmental department, KAIST, Daejeon, Republic of Korea,
Ph.D. Candidate, Civil and Environmental department, KAIST, Daejeon, Republic of Korea,
3
Associate Professor, Civil and Environmental department, KAIST, Daejeon, Republic of Korea,
Correspond to Prof. J.W. HONG (j.hong@kaist.ac.kr; jwhong@alum.mit.edu)
1

2

Keywords: SC wall, High velocity impact, Finite element analysis, Failure mode
1. INTRODUCTION
The SC walls are composed of several components suchas infill concrete slab, steel plates, and headed steel studs
that anchor the steel plate to the concrete slab. Those
structures have many advantages compared withconventional RC structure in terms of construction
efficiency and performance for the safety [1], [2], [3]. For
example, the SC walls have much higher shear strengthand resistance than the RC structure [4]. Many researchers
still elaborate in order to investigate the performance ofSC walls in different kinds of structures including nuclear
facilities as well as conventional structures [4].

solid elements while the steel studs are modeled withbeam elements. In order to describe the large deformationof steel and concrete materials, the plastic kinematic
material model and the CSCM concrete model areemployed. Contact conditions should be taken into
consideration in order to describe the interaction between
projectiles and targets as well as the contact between thecomponents of SC walls. The deformation, penetration,and perforation of the SC plates are studied by changing
the initial projectile velocity, and the critical velocity isinvestigated parametrically.
3. CONCLUSIONSThe simulation model is constructed, and the critical
impact velocity is calculated by increasing the velocity of
the impact projectile. The result is compared with thecriterion derived from the available analytical equation [5].
ACKOWLEDGEMENT
This work is financially supported by Korea Minister of
Ministry of Land, Infrastructure and Transport (MOLIT)
as「U-City Master and Doctor Course Grant Program」.

Figure 1. Schematic of SC wall.
2. NUMERICAL SIMULATION

(a)

(b)

(c)
(d)
Figure 2. Simulation models of (a) projectile, (b)concrete slab, (c) headed steel studs, and (d) steel plates.

REFERENCES[1] S. Akita, M. Ozaki, N. Niwa, I. Matsuo, and K. Hara,“Study on steel plate reinforced concrete bearing wall for
nuclear power plants (part# 2). analytical method toevaluate response of SC walls,” 16th Int. Conf. Struct.
Mech. React. Technol. (2001).
[2] M. Takeuchi, M. Narikawa, I. Matsuo, K. Hara, and S.
Usami, “Study on a concrete filled structure for nuclear
power plants,” Nucl. Eng. Des. 179 (1998) 209-223.
[3] P. F. Adams, T. J. E. Zimmerman, and J. G. MacGregor,"Design and Behavior of Composite Ice-Resisting Walls,"
ACI Spec. Publ. 109 (1987).
[4] S. Epackachi, N. H. Nguyen, E. G. Kurt, A. S.
Whittaker, and A. H. Varma, “In-Plane Seismic Behavior
of Rectangular Steel-Plate Composite Wall Piers,” J.
Struct. Eng. 1 (2014) 1-9.
[5] J. C. Bruhl, A. H. Varma, and W. H. Johnson, “Designof composite SC walls to prevent perforation from missile
impact,” Int. J. Impact Eng. 75 (2014) 75-87.

High speed impact analyses are performed using acommercial finite element code LS-Dyna, and the plastic
behavior due to the large deformation and strain rate effectare considered. As shown in Figure 2, the projectile, theconcrete slab, and steel plates are constructed using 3D― 57 ―

INFORMATION SYSTEM "COCOA" FOR EFFECTIVE EVACUATION OPERATIONBY PERSONAL IDENTITY NUMBER
Satoshi TAKATSU1, Muneyoshi NUMADA2,
Graduate Student, the University of Tokyo, Japan,
Lecturer, ICUS, Institute of Industrial Science, the University of Tokyo, Japan,
Correspond to Satoshi TAKATSU (s-taka@iis.u-tokyo.ac.jp)
1

2

Keywords: information sharing, evacuation facility, evacuee, headquarter, development of IT system, personal identity
number
1. INTRODUCTION
Information sharing among disaster headquarters and sitesof evacuation facilities is important to operate effective
evacuation facility management. But, it was difficult toobtain and share the information during earthquake
disasters. It is difficult to get even the number of evacuees
in each evacuation facility under the condition of aconcentration of a lot of evacuees.
One of the possible reason is that no effective information
sharing tools is developed. Although some IT system havebeen developed, most systems cannot manage thecondition of evacuation facilities (For example [1]). And,the developed IT system by the municipality is difficult tomaintain it without the necessary budget.
It is necessary to achieve effective and smooth operationsof the evacuation facilities by sharing the information
environments among headquarters and evacuation
facilities for the national standard system.
2. COCOA
In this study, we have developed the IT evacuation facility
system “COCOA” that can manage and support inaccordance with conditions of each evacuation facilities
by sharing the information such as the needs and demands,
the number of evacuees and those personal information.
COCOA can serve and collect the various personal
conditions and manage the information of all evacuation
facilities.
The personal identity number service (PIN service) will
start in Japan. COCOA can be used in many cases (during
staying the evacuation facility, moving to the temporary
house, condition of received administrative services, etc.)
by applying PIN service (Figure 1). For the next hugedisaster, COCOA can provide the environments forsharing the information among the evacuation facilities
and support appropriate decision-making for the person in
charge of the evacuation facility (Figure 2). Figure 3
shows the age number matrix of evacuee for each facilityon the web browser.
COCOA can manage the schedule of assigned staffs andcalculate optimum assignment schedule. The volunteer isnecessary by optimum assignment and manage theschedule. COCOA can also manage the volunteer staffs
not only government staffs.

Input:

Personal number

PeopleDB

Output:Personal info. (Name, Sex, Age, etc)

Record:Evacuation facility, timestamp

Figure 1: Information structure in COCOA
School A

School B

Headquarter

Public Building D

Public Building C

Figure 2: Information structure in COCOA

Figure 3: Age number matrix of evacuee for each facility
3. CONCLUSIONS
It is necessary to achieve smooth information sharing at
the time of the disaster. But obtaining and sharing
information was difficult during the past disasters. And
the lessons learn from the past disasters are also difficult
to update or apply to next disaster by a system. But,
COCOA is effective information tools for the evacuation
management as a Japan standard system. Japan haveexperienced for many kinds of disasters, the experiences
can be applied to other places by COCOA.
REFERENCES[1] Akari HAMANURA et al: Development of
Evacuation Support System by Considering Movement of
the User and Continuous Use in Daily Environment,
Proceedings of multimedia, dispersion, cooperation andmobile symposium, Information Processing Society ofJapan, pp.1930-1937, 2013.

― 58 ―

A ROAD TRAFFIC NOISE SIMULATION METHOD BASED ON ACOUSTIC WAVE THEORYUSING TIME DOMAIN FM-BEM
R. OKAMURA1, H. YOSHIKAWA2, T. TAKAHASHI3, and K. KASHIYAMA4
Graduate Student, Civil and Environmental, Chuo university, Tokyo, Japan,
Senior Lecturer, Applied Analysis and Complex Dynamical Systems, Kyoto university, Kyoto, Japan,
3
Lecturer, Mechanical Science and Engineering, Nagoya university, Aichi, Japan,
4
Professor, Civil and Environmental, Chuo university, Tokyo, Japan,
Correspond to Prof. K.KASHIYAMA(kaz@civil.chuo-u.ac.jp)
1

2

Keywords: boundary element method, acoustic wave theory, fast multipole method, road traffic noise, auralization
1. INTRODUCTION
As road traffic noise problems are to be serious social
problem in urban area, it is very important to evaluate
traffic noise for the planning and designing of highway.
The numerical simulation for traffic noise can beclassified into two approaches; wave acoustic theory andgeometrical acoustic theory.
This paper presents a simulation method for road traffic
noise based on acoustic wave theory. The wave scattering
problem for road traffic noise is solved by the timedomain fast multipole boundary element method (FMBEM) [1].GThe auralization is realized with the computed
sound pressure.

Figure 1. The computational model.

2. GOVERNINGE EQUATION AND NUMERICALMETHOD
The unsteady wave equation is employed for thegoverning equation. The governing equation is discretizedby the time-domain interpolation based FM-BEM
proposed by Takahashi [1], which is based on an
approximation of the integral kernels of the boundary
integral equation with interpolation functions. For the
spatial discretization, the piecewise-constant interpolation
based on triangular element is employed. On the other
hand, the piecewise-linear interpolation is employed forthe temporal discretization. The computational results are
obtained by the convolution using the impulse response
and the actual road traffic noise.
3. NUMERICAL EXAMPLE

Figure 2. Sound pressure at observation point.

To comfirm the validity of road traffic noise analysis
using the FMBEM, we calculate the three-dimensional
acoustic problem with the sound insulation wall as shownin Figure 1. The number of the boundary elements is250840, the time increment Δt is assumed to be 0.0585 ms.
The input data at the sound source point is given by theactual road traffic noise which is obtained from theobservation at the driving test course of Ministry of Land,
Infrastructure and Transport.
Figure 2 shows the comparison of sound pressure at theobservation point with or without the sound insulation
wall. The attenuation of the comparison of sound pressure
by the sound insulation wall can be clearly shown in
Figure 2. The auralization is also performed using thecomputed sound pressure. We can understand the effect ofsound insulation wall by the change of sound pressure andtone.

4. CONCLUSIONSIn this paper, we have developed a numerical simulation
method based on FM-BEM for the three-dimensionaltransient acoustic field analysis. The attenuation of soundpressure with the sound insulation wall can be verified.
The application of the present method to the actual roadtraffic problem is left for the future work.
REFERENCESG[1] T Takahashi, An interpolation-based fast-multipole
accelerated boundary integral equation method for thethree-dimensionalwave
equationୈJournalofComputational Physics Vol.258(2014), pp.809–832୊
[2] Yoshikawa H. Matsuura R. Time domainBIEM with
Lubich CQM considering arrival time of influence waves,
JASCOME, Vol.12 (2012), pp.73–78, (in Japanese).

― 59 ―

A METHOD OF MOLECULAR COMPUTATION OF PHOTO-RESPONSIVE POLYMERREINFORCED WITH NANOPARTICLESJOONMYUNG CHOI1, JUNGHWAN MOON1 and MAENGHYO CHO2
1

Ph. D. Student, School of Mechanical and Aerospace Engineering, Seoul National University, Seoul, South Korea
2
Professor, School of Mechanical and Aerospace Engineering, Seoul National University, Seoul, South Korea
Correspond to Prof. MAENGHYO CHO (mhcho@snu.ac.kr)

Keywords: Molecular dynamics simulation, Photo-responsive polymer, Nanocomposites, Nano-mechanical modeling,
Gold nanoparticle, Carbon nanotube
1. INTRODUCTION
In this paper, we introduce the molecular-levelcomputational modeling method, which for the optomechanical design and analysis of photo-responsive
polymer (PRP) and polymer nanocomposites. PRPs havetheir unique characteristics in that the material reversibly
deforms up to centimeter-scale by irradiation of specificwavelength of light [1]. The main cause of such photoresponsibility of the material is that the photo-reactive
mesogen molecules inside of the polymer networks
directly modulate the alignment of the polymer
microstructures under light/thermal exposure. In thepresent work the effective photostrain and thecorresponding thermal and mechanical stability of thematerial are quantitatively derived according to thephoto-isomerization ratio of the reactive mesogen
molecules [2], by adapting photo-reactive potential to theclassical molecular dynamics (MD) simulations.
Together with the molecular computation model forbulk responsibilities of PRP structure, the present work
also covers the improvement of the thermal andmechanical properties of the material by inserting nanosized gold particles or carbon nanotubes. Due to thestrong filler-matrix interaction and its extremely highsurface-to-volume ratio [3], one can expect theextraordinary properties as well as multi-functionalitiesof the nano-sized filler embedded composite system.
2. COMPUTATION MODELS
To simulate the photo-responsibility of PRP structure
and its nanocomposites, azobenzene-contained liquid
crystalline polymer molecules are generated in silico andallocate along nematic direction. Regarding thenanocomposites system, silicon carbide (zerodimensional) and single-walled carbon nanotube (onedimensional) filler are positioned in the center of the unit
cell as shown in Figure 1. The polymer matrix aretightly crosslinked and thermodynamically equilibratedusing MD. The photoactive potential is selectively
applied in the unit cell to estimate the photostrain andother thermo-mechanical properties of the material
according to the photo-isomerization ratio of azobenzene
molecules. The elastic stiffness and the thermal stabilityof the considered unit cells are derived through uniaxial
tensile loading simulation under NVT and heating-up
simulation under NPT ensembles, respectively.

Figure 1. Equilibrated molecular configurations of the
PRP/nanoparticle composite structures constructed in the
MD simulation (left) and their modeling process (right).
3. CONCLUSIONS
It is verified that the considered 0-D and 1-D
nanocomposite models show significantly enhanced
mechanical and thermal properties without loss of photodeformation characteristics of polymer domain. The
effect of interfacial region between the nanofiller and thesurrounding matrix and the nanoparticle particle size
effect also be characterized.
This work was supported by the National Research
Foundation of Korea (NRF) grant funded by the Korea
government (MSIP) (No. 2012R1A3A2048841)
REFERENCES[1] Y. Yu, M. Nakano, T. Ikeda, Nature 425 (2003) 145.
[2] J. Choi, H. Chung, J.-H. Yun, M. Cho, Appl. Phys.
Lett. 105 (2014) 221906.
[3] J. Choi, H. Shin, S. Yang, M. Cho, Compos. Struct.
119 (2015) 365-376.

― 60 ―

STABILIZED FINITE ELEMENT METHOD FOR WIND FLOW ANALYSIS WITH PLANT
CANOPY
S. HORIIKE1, and K. KASHIYAMA2
1

Graduate student, Department of Civil and Environmental Engineering, Chuo University, Tokyo, Japan
2
Professor,Department of Civil and Environmental Engineering, Chuo University, Tokyo, Japan
Correspond to Prof. K. KASHIYAMA(kaz@civil.chuo-u.ac.jp)

Keywords : LES, plant canopy model
1. INTRODUCTION
In recent years,G the heat-island phenomenon is
becoming serious problems in the urban area by theeffect of rapid urbanization. The major cause of the heatisland phenomenon is based on the decrease ofvegetation, and thermal storage effect by alteration ofground and structures. In order to estimate the extent of
the environmental issues quantitatively, it is develop an
accurate numerical method. The finite element method issuitable method for the thermal environment in urban
area, since the finite element method can easily treat thegeography and buildings.
This paper presents a numerical method for theevaluation of environmental flow with plant canopy by astabilized finite element method. The one-equation
model of LES [1] is employed for the turbulence model.
2. NUMERICAL METHODSIn this study, we deal with the compressible viscous
fluid. The incompressible Navier-Stokes equation andequation of continuity are employed for the basic
equation of flow field. For the turbulent model, the oneequation model based on LES [1] is employed. In order
to consider the effect of plant canopy, the plant canopy
model [2] is employed.
For the discretization in space, we employed thestabilized finite element method based on the
SUPG/PSPG method. The Crank -Nicolson method is
employed for the discretization in time. The parallel
computing based on MPI is performed. The element-byelement Bi-CGSTAB2 method is employed for solving
the discretized simultaneous equations.
3. NUMERICAL EXAMPLE
The present method is applied to a numerical example
containing a plant canopy to investigate the validity andefficiency of this method. The computed results arecompared with existing experimental results [3].
Figure 1. shows the computed domain.G The inlet
velocity is assumed to the values obtained from theexperimental data [3], and the outlet velocity is free
outflow boundary condition. The upper and side is slip
boundary condition. The parameter for plant canopy
model are assumed as; the drag coefficient of tree crown
C f 0.8 , the leaf area density a 17.93 [m2/m3] . The
minimum mesh size and the time increment are assumed
to be 'x

0.06 and 't 1.0 u103 . Figure 2. shows

Figure 1. Numerical ExampleG

Figure 2. Velocity distribution(y=0.6[m], y=0.8[m])the velocity distribution at y=0.6, y=0.8. The computed
result obtained by the present method is good agreementwith the experimental value.
3. CONCLUSIONSThis paper presents a numerical method based on thestabilized finite element method for the wind flow
considering the plant canopy model. From the results
obtained in this paper, it can be concluded that thepresent method is an useful tool for the wind flow
analysis. The consideration of the radiation damping by
the plant canopy is left for the future work.
REFERENCES[1] Deardorf,J.W. Stratocumulus-capped mixed layers
derived from a three-dimensional model, BoundaryLayer Meteorology, (1980) 18 , pp.495-527[2] H.Hiraoka, Modeling of a one-equation LES for plant
canopy flows, J.Environ. Eng. AIJ, Vol.74 㸪 No.639,(2009)603-612. (in Japanese)[3]GM.Ohashi, A numerical simulation on air flow aroundmodel tree, Building Research Institute, Ministry ofConstruction, (1998)713-714. (in Japanese)― 61 ―

NUMERICAL INVESTIGATION OF NATURAL CONVECTION EFFECTS COUPLED WITH
THERMOACOUSTIC OSCILLATING FLOW IN DIFFERENT GRAVITY INCLINATION
J.Q. HAN1, J.G. CHEN2, H.Y. YUAN2, and J.P. ZHANG1
1

Post doctor, Institute of Public Safety Research, Tsinghua University, Beijing, China,
Professor, Institute of Public Safety Research, Tsinghua University, Beijing, China,
Correspond to Prof. J.G. CHEN (chenjianguo@mail.tsinghua.edu.cn)

2

Keywords: thermoacoustic, natural convection, gravity effects, oscillating flow
1. INTRODUCTION
Pulse tube refrigerator has many advantages overconventional refrigerators, such as no moving part, lowcost, high reliability and less mechanical vibration.
Therefore, it is recognized as one of the most promisingrefrigerators in the future. There exists gravity property inpulse tube refrigerators, as shown in figure 1, whichincarnates that the gravity affects thermoacoustic engines
efficiency when the hot side of the pulse tube is placed in
the lower than the cold side. This problem attracts many
researchers, including Fujimoto, Johnson, Thummes,
Kasthurirengan and others. Previous research results arealmost experiments. Quantitative results and mechanism
knowledge which need theoretical analysis and numericalsimulation are seldom. There are many factors that affectthe efficiency pulse tube refrigerator, including the gravity
angle, the oscillating flow frequency, the temperature
difference between the hot and cold side, the pressure
amplitude of oscillating flow, and so on. This article
researches the coupling mechanism between oscillating
flow and natural convection in pulse tube with the changeof angle by numerical simulation.

Figure 1. Thermoacoustic pulse tube refrigerator
2. RESULTS AND DISCUSSIONS
We study numerically oscillating flow coupled withnatural convection in different inclination angle by thefinite volume method with the SIMPLE algorithm.
Boundary layer mesh is used to get accurate details of theflow field. α is the angle between gravity inclination andtemperature gradient. �� is the temperature difference
between cold side and hot side. Nu number characterizes
convective heat transfer intensity. The development of
natural convection is governed by the Navier-Stokes
equations for compressible fluid.
Figure 2 is the comparison of oscillating flow Nu andnatural convection Nu in different gravity angles at �� �
����. Nunatural is the Nu number for natural convection
and Nuoscillating is the Nu number for oscillating flow
coupled with natural convection. In these two cases, thetemperature difference, geometry of structure and otherconditions are same.

Figure 2. Oscillating flow Nu and natural convection Nu
in different gravity angles at �� � ����.

3. CONCLUSIONS
We find there is the intersection of two Nu curves,which is corresponding to a critical angle. When it is lessthan the critical angle, the oscillating flow Nu is smaller
than natural convection Nu and oscillating flow weakens
natural convection; and when it is greater than the criticalangle, oscillating flow Nu is larger than natural convection
Nu, and oscillating flow is enhanced by natural convection.
In this article, for the first time, the couplingmechanism of the oscillating flow and the naturalconvection with the change of inclination angle isrevealed. For the latter part of the study, this part work
provides a reference basis.
4. REFERENCES[1] S.Fujimoto, Y.M. Kang, and Y. Matsubara,
Cryocoolers, 109(1999) 213-220[2] S. Kasthurirengan, Advances in Cryogenic
Engineering, 9a (2004) 1474–1481.
[3] S. Kasthurirengan, G. Thummes, and P.T. Zubkov,
Advances in Cryogenic Engineering, 1 (2000)136-140.
[4] M. Shiraishi, Cryogenics, 44 (2004) 101-107.
[5] M. Shiraishi, Advances in Cryogenic Engineering, 49a(2004) 1560–1567.
[6] G. Thummes, Cryocoolers, 9 (1997) 393–402.
[7] L.M. Zhang, J.Y. Hu, E.C. Luo, Cryogenics, 51 (2011)85–89.

― 62 ―

CODE2015

Technical Program – Day 2

Session 6

STRUCTURAL TOPOLOGY OPTIMIZATION OF ELECTROMAGNETICVIBRATION ENERGY HARVESTERS
J. LEE1, and S.W. Yoon2
1

Assistant Professor, School of Aerospace and Mechanical Enginnering, Korea Aerospace University, Gyeonggi-do,
South Korea
2

Assistant Professor, Department of Automotive Engineering, Hanyang University, Seoul, South Korea
Correspond to Prof. J. LEE (wookslee@kau.ac.kr)

Keywords: Topology optimization, Energy harvesters, Permanent Magnet
1. INTRODUCTION
This work presents topology optimization ofelectromagnetic vibration energy harvesters [1]-[2]. The
harvester transduces the mechanical vibration into electricenergy through the magnetic field based on theelectromagnetic induction. The harvester is composed ofa mechanical coil spring, wires and permanent magnets(PMs) with or without a magnetic back iron. In this work,
the configuration and geometries of the permanent magnet
with or without back irons are designed using the topology
optimization approach. For both problems with andwithout back irons, the magnetization direction of thepermanent magnet is also designed. The optimization
objective is the maximization of the root-mean-square(RMS) value of the harvested output voltage, which will
determine the amount of the harvested energy. The output
voltage is calculated using the Lagrange polynomialinterpolation of the flux linkages, which can be calculated
from the magnetic field distribution. For calculating themagnetic field distribution, the Maxwell’s equation issolved using the finite element analysis for the given
mechanical vibration. The sensitivity of the objectivefunction is derived using the adjoint variable method. The
volume of the permanent magnet and back irons areconstrained as the given values. For the topology
optimization the conventional nodal-density based
topology optimization is applied using the commercial
software COMSOL.
2. OPTIMIZATION RESULTThe proposed topology optimization approachsuccessfully finds the optimal design of the harvester withor without back iron structure. The haversters are modeled
as axi-symmetric two dimentional one. The design resultwith and without back iron is obtained. Fig. 1 shows thedesign result with back iron when the volume ofpermanent magnet is 25% and the back iron volume is also
25% of the whole design domain. The optimal design
result is compared with the initial benchmark design. The
comparison results shows the huge increase of the output
voltage. This enhancement comes from the optimal
location and shape of the harvester components, which
maximize the time variation of the flux linkages. In

addition, the effect of the volume constraint for thepermanent magnet and/or back iron is also investigated.

Fig. 1. Topology optimization result of electromagnetic
vibration energy harvesters.
3. CONCLUSIONSThis work presents the topology optimization ofelectromagnetic vibration energy harvesters. The
permanent magnet with or without back iron is
successfully optimization to maximize the harvester
performances. In the future, the optimal design will befabricated, and the performance enhancement will bevalidated by the experiment.
REFERENCES[1] J. Lee, S. W. Yoon, Topology design optimization ofelectromagnetic vibration energy harvester to maximize
output power, Journal of Magnetics, 18(3) (2013) 1-6.
[2] J. Lee, S. W. Yoon, Optimization of Magnet and BackIron Topologies in Electromagnetic Vibration Energy
Harvesters, IEEE Transactions on Magnetics, 51(6)(2015).

― 63 ―OIL SPILL SPREAD SIMULATION MODEL ON COMPLEX TOPOGRAPHY BASED ON
HYDRODYNAMICS
Boni Su1, Hong Huang2, and Yuntao Li3
Ph. D. Candidate, Institute of Public Safety Research, Tsinghua University, Beijing, China,
2
Professor, Institute of Public Safety Research, Tsinghua University, Beijing, China,
3
Lecturer, School of Mechanical and Transportation Engineering, China University of Petroleum, Beijing, China,
Correspond to Prof. H. Huang (hhong@tsinghua.edu.cn)
1

Keywords: Oil spill, Simulation Model, Complex Topography, Hydrodynamics
1. INTRODUCTION
Accidents of oil leakage occur frequently duringtransport and storage. Such accidents waste oil resources,
pollute environment, and may cause serious fire.
Numerical simulation of oil spill spread is helpful for risk
analysis and accident prevention.
Most of studies on oil spill spread paid attention towater surface or flat ground [1]. However, experiments
showed that ground slope influence the spreadsignificantly [2], indicating that topography is an
important factor in fuel spread. Therefore, it is almost
impossible to simulate oil spill spread accurately underreal accident scenarios due to the lack of simulation modelon complex topography.
In this paper, an oil spill spread simulation model oncomplex topography is established based onhydrodynamics. The model is applied to an actual oil tankarea for oil leakage simulation.
2. METHODS
SWEs (Shallow Water Equations) are widely used forflood simulation in hydraulics. In this paper, oil spill
spread simulation model is established based on adiffusive wave approximation model, LISFLOOD-FP [3].
Many details were improved to make the model suitablefor oil spill spread simulation. Firstly, the calculation
formula of flow velocity is modified from Manning
Formula to Darcy-Weisbach Equation. Reynolds number
in oil spill spread is usually small (due to high viscosityand low velocity) and maybe not in turbulence regime.
Darcy-Weisbach Equation is suitable for all flow regimes,
while Manning Formula is only applicable in complete
turbulence regime. In this paper, Darcy friction factor wasdetermined based on formulas in reference [4]. Secondly,
additional limit condition at spread border is introduced toobtain correct spread area. Flow is allowed to flow out agrid cell only when depth is greater than a threshold [5].
3. RESULTS AND DISCUSSIONSAn oil tank area is selected as the study area. Figure 1(a)shows the region of study area, and Figure 1(b) shows the
GIS (Geographic Information System) data.
An accident scenario is simulated. In this scenario, oiltanks #9, #10, #11, and #12 leak at the same time due toexternal damage. The tanks contain crude oil. Figure 2
shows the simulation results, which are the spatialtemporal distributions of oil depth.

Oil gradually filled the pool, which is designed to stop
fuel spread. Due to the limit of pool volume, oiloverflowed form the pool and started spread to low-lying
areas. The spread area covered most of the area, large
scale fire may occur if oil is ignited, which is very
dangerous.

(a) Region
Figure 1. Study area(b) GIS Data

Figure 2. Simulation results
4. CONCLUSIONSIn this paper, an oil spill spread simulation model oncomplex topography is established based onhydrodynamics. Many details were improved to make themodel suitable for oil spill spread simulation. The model
is applied to an actual oil tank area for oil leakage
simulation, and the spatial-temporal distributions of oil
depth were obtained. The simulation model is helpful forrisk analysis of overland oil spill.
REFERENCES[1] D.P. Hoult, W. Suchon, The spread of oil in a channel,
Massachusetts Institute of Technology, Cambridge, 1970.
[2] C.S. Simmons, J.M. Keller, J.L. Hylden, Spills on Flat
Inclined Pavements, Pacific Northwest National
Laboratory, Richland, 2004.
[3] P.D. Bates, A.P.J. De Roo, A simple raster-based model
for flood inundation simulation, J. Hydrol. 236(2000) 5477.
[4] B.C. Yen, Open Channel Flow Resistance, J. Hydraul.
Eng. 128(2002) 20-39.
[5] M. Hussein, M. Jin, J.W. Weaver, Development andverification of a screening model for surface spreading of
petroleum, J. Contam. Hydrol. 57(2002) 281-302.

― 64 ―

A FINITE ELEMENT STUDY ABOUT MECHANICAL FRUSTRATION OF NEMATIC SOLIDS
INDUCED BY DISCLINATION DEFECTHAYOUNG CHUNG1, JAESUNG PARK1, JOONMYUNG CHOI1, JUNG-HOON YUN1,
and MAENGHYO CHO2
1

Ph.D. student, Mechanical engineering, Seoul National University, Seoul, South Korea,
2
Professor, Mechanical engineering, Seoul National University, Seoul, South Korea,
Correspond to Prof. M. Cho (mhcho@snu.ac.kr)

Keywords: Nematic solid, Disclination defect, Finite element, Geometric nonlinearity
1. INTRODUCTION
A nematic solid, wherein a uniaxial long-range order isincorporated into a polymer via crosslinking of rigid rod
and flexible chain, has been envisaged as a multi-physicalsmart material as macroscopic shape change is induced byeither mechanical load or order-changing stimuli such as
light and temperature [1]; upon exposure to stimuli thatchange alignment condition of liquid crystalline
chromophores, the structure either shrinks/expands
following the principal direction and the of the shape ofthe microstate. Nematic director n, which terms an
ensemble average of the rigid rods, is thus determines thedeformed shape of the structure. For instance, when the
director is deviated from the longitudinal direction of thecantilevered thin strip, the structure twists, and even bend
away from light.
Accordingly, nontrivial distribution of the director overthe structure generates corresponding stress field that
oftens induces mechanical frustration. A disclination
defect, by which a arbitrary Schileren texture found in theliquid crystal is parameterized, is therefore pivotal tounderstand order-induced deformation.
In this work, we numerically investigate thedeformation and the mechanical frustration of the nematicsolid generated by disclination-like texture. A core
position and the strength of the defect is blueprinted
arbitrarly; a corresponding stress field generated byphotomechanical constitutive relation and its effect on the
structure is analyzed by facilitating nonlinear finite
element
methodincorporating
photomechancal
discussions [2].
2. METHODS
A nematic solids with liquid crystal chromophore
deforms as the uniaxial order of rigid constituents either
increases or decreases. Such mechanical phase behaviorof projecting microscopic conformation to macroscopic
structure is well explained by reference [2]The constitutive equation is incorporated into
geometric nonlinearity with corotational consideration.
Element Independent Corotational (EICR) formulation,
which utilizes displaced element to extract in/out rigid
body displacement and rotation. Discrete Kirchhoff
Triangle (DKT) and Optimal membrane triangle (OPT)elements are used to construct linear 18-dof element.
Distribution of the nematic director n(x) is assumed toknown a priori and be blueprinted before the change oforder. Some examples of disclinations are tested in thepresent work; the positions of the core and strength arealtered arbitrarily, where the distributions at the vicinal of
the core are computed by 1st order Fredericks distortion
energy whenever required.
3. RESULTS
A light-induced deformation topographies of circular
nematic solids with disclination strength m and center core
position are illustrated in Fig. 1. Monopole (m=1) and
quadpole (m=-1) geometry are exhibited, and both agree
well with the experimental results. Geometric nonlinearity
is found to be crucial to correctly figure out such topology
evolution, as linear counterpart does not relax mechanicalfrustration imposed by disclination, and thus possesses
overestimated membrane strain energy.

Figure 1. Deformed pattern of photo-responsive nematic
liquid crystal with texture marked by m.
3. CONCLUSION
In this work, we simulate deformation of nematic solid
induced by mechanical frustration with blueprinted
disclination; by facilitating finite element formulation thatcombines nematic constitutive relation with geometric
nonlinearity, effects from arbitrary texture are fairly
simulated when compared to experiments.
ACKNOWLEDGEMENT
This work was supported by the National Research
Foundation of Korea(NRF) grant funded by the Korea
government(MSIP) (No. 2012R1A3A2048841)
REFERENCES[1] M. Warner, E. M. Terentjev, Liquid Crystal Elastomers,
second ed., Clarendon Press, 2006.
[2] H. Chung, J. Choi, J.-H. Yun, M. Cho, Phys. Rev. E(2015) 042503― 65 ―

PEDESTRIAN SPEED CHANGE MODEL FOR SAFETY ASSESSMENT AT SIGNALIZEDCROSSWALKS
M. IRYO1, W.K.M. ALHAJYASEEN2
Associate Professor, Institute of Industrial Science, the University of Tokyo, Tokyo, Japan,
Assistant Professor, Department of Civil Engineering, King Fahd University of Petroleum and Minerals, Dhahran,Saudi Arabia,
Correspond to Assoc. Prof. M. IRYO (m-iryo@iis.u-tokyo.ac.jp)
1

2

Keywords: Traffic Safety, Pedestrian Behavior, Empirical Data Analysis, Traffic Simulation

2. PEDESTRIAN SPEED CHANGE EVENTS AS
CRITICAL MANEUVERS
Pedestrians’ sudden speed change is an important factor
for safety assessment. Figure 1 shows an example ofpedestrian trajectory. If there is a turning vehicle
approaching to the conflict area, the driver may behave
based on the expected arrival time of pedestrian at the area.
As he/she suddenly changed the speed in the middle ofcrosswalk, the difference between the expected and actual
pedestrian arrival time, Tdif, becomes approx. 2 seconds,
which is significantly large in terms of conflict avoidance.
Empirical observation data show that pedestrians tend
to change their speed when they are approaching theconflict area and when there is a large gap between current
speed of the subject pedestrian and the necessary speed tocomplete the crossing. Further in-depth analysis isconducted to quantitatively explain the probability of
speed change considering location, timing and otherrelated variables.
3. SIMULATION INTEGRATION FOR SAFETYASSESSMENT
The developed model is applied in a simulation modelthat contains a comprehensive turning vehicle maneuvermodel [3]. Pedestrian-vehicle conflicts are generated by
the traffic simulation.
The traffic simulation model can derive surrogate safetymeasures such as Post Encroachment Time and vehicle
speed at a conflict point, in order to indirectly evaluate theprobability and the severity of pedestrian-vehicle conflicts.
Sensitivity of the pedestrian maneuver upon surrogate
safety measure is analyzed through the simulation.

Far-side
Pedestrian

35
Spacing from near-sidecrosswalk
Spacing(m) edge

1. INTRODUCTION
Pedestrians are vulnerable road users. In Japan, morethan one-third of the fatalities in traffic crashes arepedestrians. To evaluate alternative countermeasureseffectively, recently traffic simulation is considered as oneof the powerful decision support tools. So far, most of theexisting simulations assume that pedestrians walk at aconstant speed and follow the traffic rules when they crossroads. In reality, the pedestrians may suddenly startrushing in order to complete the crossing or for otherreasons. These behaviors are critical for collision
avoidance.
The objective of this study is to integrate a detailed
pedestrian maneuver model at signalized crosswalks into
the traffic simulation model for assessing traffic safety.

30
Expected pedestrianarrival time25
Actual pedestrianarrival time

20
15
10

Tdif

5
0

Near-side

Pedestrianvehicle
conflict area

-5

30

35

40

45

50

Elapsed time from onset of pedestrian green(second)Turning vehicle

Figure 1. Example of a pedestrian speed change event anda possible pedestrian-vehicle conflict
4. CONCLUSIONSThis presentation demonstrated the impact ofpedestrian speed change maneuver on traffic safety. The
proposed model can be used for safer crosswalk/traffic
signal design and the development of vehicle safetyassistance systems.
ACKNOWLEDGEMENTS
This research is supported by the Deanship of Scientific
Research (DSR) at King Fahd University of Petroleum &Minerals (KFUPM) through project No. IN141020 as well
as JSPS KAKENHI Grant No. 15H05534. The authors are
grateful to Prof. Hideki Nakamura for providing video
survey data.
REFERENCES[1] U. Shahdah et al., Application of traffic
microsimulation for evaluating safety performance ofurban signalized intersections, Transp. Res. C, 60 (2015)
96-104.
[2] M. Iryo-Asano et al., Analysis and Modeling ofPedestrian Crossing Behavior during Pedestrian Flashing
Green Interval, IEEE Trans. on ITS, 16(2) (2015) 958-969.
[3] M.T. Dang et al., Development of a Microscopic
Traffic Simulation Model for Safety Assessment at
Signalized Intersections, Transp. Res. Rec., 2316 (2012)
122-131.

― 66 ―

NON-INTRUSIVE REDUCTION METHOD WITH CONNECTED ELEMENTS
E. KIM1, and M. CHO2
1

Ph. D. candidate, Mechanical and Aerospace Engineering, Seoul National University, Seoul, Korea,
2
Professor, Mechanical and Aerospace Engineering, Seoul National University, Seoul, Korea,
Correspond to Prof. M. CHO (mhcho@snu.ac.kr)

Keywords: Maximum Reduced Order Model, Structural dynamics, Optimization, Nonlinear transient analysis,
Geometrical nonlinearities, Multibody Dynamics, Absolute Nodal Coordinate Formulation
1. NON-INTRUSIVE
REDUCED ORDER MODELING
Various reduced order modeling techniques have beenstudied to solve the geometrically non-linear transient
response analysis. The nonlinearities cause the iterativeevaluations of the internal force, which makes the
reduction process difficult. If we reduce the internal force
itself, the equivalence of the system suffers. If we reduce
the system without the internal force reduction, theefficiency does not rise much.
As shown in the eq.(1), a non-intrusive reduction
modelling technique uses the polynomial formulation todescribe non-linear internal force terms. With the equation,
the construction time of the reduced model increases withcubic order of the system size.[1,2]M ij q j + Cij q j + Gi (q) = Fi
(3)
Gi (q) = K ij(1) q j + K ijl(2) q j ql + K ijlp
q j ql q p

This reduced order model can be effectively included in
the optimization process. See Fig. 1.
2. NUMERICAL EXAMPLE
A cantilever beam is subjected to the dynamic loadingat the tip. A 9-node nonlinear shell element with 6-DOF
per node is used. The length and width of the beam are
0.5m and 0.05m, respectively. The number of elements isvaried from 8(270 DOFs) to 48(1,350 DOFs). The NICE
reduced model is established and the transient analysis isconducted under the perturbed loading condition. Fig. 2.
shows the effectiveness of the proposed method. The
reduced model is 29 times faster for 8-elements case and37 times faster for 48-elements case.

(1)

In this work, a Non-Intrusive reduction method with
Connected Element (NICE reduction method) is proposed.
The element connectivity is used to construct the nonintrusive models. Then, the system equation is reduced toeq.(2) with proper orthogonal modes F .

FT M F a +FT C F a +FT G(q) = FT F
qi = Fij a j(2)
Figure 2. The computation time for the transient analysis
3. CONCLUSIONS
A new reduced modelling technique (NICE reduction)
for the geometrically nonlinear transient analysis is
proposed. Non-intrusive modelling technique wasmodified using the information on the element
connectivity. The NICE reduction method is useful for thedesign optimization of the nonlinear transient system. The
larger the system, the more effective the method becomes.
AKNOWLEDGEMENT
This work was supported by the National Research
Foundation of Korea(NRF) grant funded by the Korea
government(MSIP) (No. 2012R1A3A2048841).

Figure 1. Flow chart for the optimization of transient
dynamic system based on NICE reduction method.

REFERENCES[1] M.P. Mignolet, A. Przekop, S.a. Rizzi, S.M.
Spottswood, J. Sound Vibration, 332 (2013) 2437-2460.
[2] F. Adrian, L.D. Tran, H.G. Matthies, R. Ohayon,
Comput. Mech., 55 (2015) 327-344.

― 67 ―

STRESS CONSTRAINTS FOR DUCTILE AND BRITTLE MATERIALS IN FLUID-STRUCTURE
INTERACTION SYSTEM
G.H. Yoon1
1

Associate Professor, School of Mechanical Engineering, Hanyang University, Seoul, Korea
Correspond to Prof. G.H. YOON (ghy@hanyang.ac.kr)

Keywords: Stress-based topology optimization, fluid-structure interaction, monolithic approach
1. INTRODUCTION
This study developed a new stress-based topology
optimization method (STOM) for a steady state fluidstructure interaction (FSI) structure minimizing thevolume, subject to the local stress constraints for brittle
and ductile materials. Despite some researches on thesubject of STOM, perplexing optimization issues related
to stress-based topology optimization for FSI
multiphysics systems still are. Serious issues whencreating a effective topology optimization for an FSI
structure include: the interpolation method between thefluid equation and the structure equation with respect tolocally defined design variables; the mutual multiphysics
coupling boundary conditions at dramatically evolving
interfacing boundaries; and a clear clarification of thegoverning equations and the interaction boundary
conditions for varying intermediate design variables asshown in Figure 1. Along with these three issues, which
are related to multiphysics equations, there are threeimportant considerations related to the STOM: the stresssingularity issue, the issues of multiple constraints and thehighly nonlinear behavior of the stress constraints.
Particularly this problem becomes serious when the
failure constraints for ductile and brittle materials are
considered. To resolve all of the aforementioned issues,
we applied a monolithic analysis, integrating the qprelaxation method and the global p-norm approach. And
we employ the differentiable max and min operators. In
our present method, we were able to create optimallayouts that minimize the volume constraining local stressvalues for a steady state fluid and structural interaction
system for brittle and ductile materials.
Before deformation

CS ? C fCSf  0approach is employed (See [1] for more details).


T

T

X

0



 F   v T F d     v v F d    p nd
T



T

X

0

T

f

p

0



0




p



(1)

*

*

  p (   v ) F d  0

(2)

T

x

0



  S  T d
T

0

s



    F  S ( u,  u )  p F d      F  u   p F d 
T

0



T

T

0



( 3)

X

All the notations are explained in the reference [1].
There are many failure theories for brittle and ductile
materials such as Tresca failure theory, Columb-Mohr
theory, Mohr theory and Drucker-Prager yield theory. In
order to consider these failure theories in topology
optimization, the differentiable max and min operators
should be developed.
3. CONCLUSIONSThe present research adopts the monolithic approach for
FSI analysis for topology optimization and develops somedifferentiable failure theories for ductile and brittle
materials. With the developed approach, the failure
constraints for various materials can be considered.
[1] G.H.,Yoon, Comput. Methods Appl. Mech. Engrg.
278 (2014) 499–523.

S  

After deformation

Structural domain  v  ( v  F  v ) F d 

Fluid domain

Figure 1. The monolithic analysis and design method fora steady-state fluid-structure interaction problem without
explicit interaction boundaries: material values forstructural and fluid domains[1].
2. MONOLITHIC APPROACH FOR FSI AND
FAILURE THEOREISIn the present approach, the following monolithic― 68 ―

CODE2015

Technical Program – Day 2

Session 7

OPTIMAL DESIGN OF TIMBER THROUGH COLUMN
H. SATO1, M. KOSHIHARA2 and T. MIYAKE3
Research Associate, IIS, the University of Tokyo, Tokyo, Japan,
2
Professor, IIS, the University of Tokyo, Tokyo, Japan,
3
President, Nihon System Sekkei Architects & Engineers Co., Tokyo, Japan
1

Keywords: Traditional Timber Construction, Traditional Townhouse, Seismic Analysis, Parametric Study
1. INTRODUCTION
This paper presents the study of analytical evaluationof timber frame with through column for the optimal
design of timber through column.
Japan has a long history of earthquakes. The oldtimber structures in Japan have suffered great damage
caused by strong earthquakes, because many of the
traditional structures often have insufficient earthquakeproof performance. Besides, many of traditional timber
townhouses in each area have same construction systemor same structural elements. If the structural performanceof their system and elements is clarified, the technique ofearthquake-proofing suitable for those buildings can be
examined. Therefor it is important to clarify structural
performance of their elements in each historical area.
2. SEISMIC ANALYSIS
To evaluate the structural performance of throughcolumns, the analyses were performed using the mass
system model added the through columns, as shown intable 1 [1]. The structural elements in this analysis were
mud walls at the second floor and through columns.
The input waves of the analysis were five simulated
earthquake motions equivalent to design earthquake
ground motion based on the Japanese Code. They weremodulated to the standard level using coefficient 0.85.

The stiffness at the second floor and the column
dimensions were used for parameters. To verify the
optimally design, the failure mode was established twocase as follows; (1) the collapse caused by few quantities
of walls and (2) the bending fracture of through column
caused by much quantities of walls.
As results, the optimal design of timber through column
considering the wall quantity at the second floor wasverified, as shown in figure 1.
When the limit of collapse was established at the 1/10
rad., the optimal wall quantity was approximately 2.5 P
or more. When the limit of collapse was established at
the 1/30 rad., which was the safe limit for new
construction, the optimal wall quantity did not exist
without through column more than 210 mm. In addition,
the column dimension and the failure mode by bending
fracture of through column have correlation.
(1)

135×135
×150
×180
×210
150×150×180×210
180×180
210×210240×240300×300(2)

1/20rad.

1/30rad.

0

Table 1. Parameter of mass system model

1/10rad.

2

4

6

8

10

12

14

16

18

20

Optimal wall quantity (P)Figure 1. Optimal wall quantity at the second floor

Model

Structuralelements
nc
Weight (tf)nw (P)
ColumnSection (mm)Species
Input level

K1 = 0 (Mud walls at 1st floor)
K2 = nw×P (Mud walls at 2nd floor)
I×nc (Through columns)
14
M1F: 2, M2F: 5, MRF: 5
1 – 20 (1 P = 0.91 m)
135×135, ×150, ×180, ×210
150×150, ×180, ×210
180×180, 240×240, 300×300
Japanese cedar
fb: 70 kgf/cm2，E: 70 tf/cm2
Japanese Code×1.0

3. CONCLUSIONSIn this model with actual situation of townhouse, theoptimal design of timber through column considering thewall quantity at the second floor was verified.
1) When the limit of collapse was established at the 1/10
rad., the optimal wall quantity at the second floor was
2.5 P or more independently of the column dimension.
2) When the limit of collapse was established at the 1/30
rad., the optimal wall quantity did not exist withthrough columns of small dimension.
3) The column dimension and the failure mode bybending fracture of through column have correlation.
REFERENCES[1] H. Sato, M. Koshihara, T. Miyake, Horizontal load
capacity of sashikamoi frame with through column Part 2
Effects of seismic performance of second floor on
through column, Summaries of Technical Papers ofAnnual Meeting Architectural Institute of Japan, C-1(2012) 461-462― 69 ―

NUMERICAL ANALYSIS OF 2D-INFILLED FRAME WITH INTERFACE MATERIALSBETWEEN REINFORCED CONCRETE FRAME AND INFILL WALL
Y. PATEL1 and K.MEGURO2
Graduate Student, Department of Civil Engineering, University of Tokyo (UTokyo), Tokyo, Japan,
Director/Professor, International Center for Urban Safety Engineering (ICUS), Institute of Industrial Science (IIS),
UTokyo , Tokyo, Japan
1

2

Keywords: Infilled Frame, Monotonic Loading, Interface Material, Applied Element Method
1. ABSTRACT
Infilled structures are widely used for construction in
many parts of the world due to social acceptability andease in availability of construction materials. Observationfrom last major earthquakes and previous studies havedemonstrated that infill wall in reinforced concrete frame
causes an increase in initial stiffness and strength.
However, it has been also observed that infill walls are thefirst element to be damaged due to brittle nature andincrease of attracted force during earthquake (increase of
attracted forces often associated with the change indynamic properties due to increase in stiffness). Infill wall
has four types of damage; crushing at the corner, diagonaltension failure, shear slide failure and out-of-plane failure.
Current study is focused on the in-plane characteristic ofinfilled frame. The current aim of the study is to propose
a simplified method for modification in infilled frame to
reduce in-plane damage in masonry wall by the use of asoft material at the side interface (column and infill wall)
and weak mortar at the top interface (beam and infill wall).
Effect of proposed method were understood on nonengineered and engineered infilled frame. And, underdisplacement controlled loading their shear resistant,storey drift and displacement ductility of non-engineered
and engineered frames were observed.
2. INTRODUCTION
During an earthquake disaster, loss of non-structuralcomponent in a building may easily surpass thereplacement cost of building [1], due to loss of business,
downtime and inventory. In addition, many constructions
in developing countries are not regularized and done
without proper monitoring, which increases the risk of life
safety as well as economic losses, especially in the nonengineered frames. There has been many research
focusing on strengthening of reinforced concrete infilled
frames but very few interest has been shown in the
reduction of damage to infill walls and harnessing that
intact infill wall strength. With the introduction ofinterface material in the infilled frame, even if reinforced
concrete structural member suffers damages, infill wall
can support the upper storey load and could effectively
mitigate the partial or full collapse.
Numerical simulation of engineered and nonengineered frames were analyzed using Applied Element
Method (AEM) [2]. AEM divides the infilled frame in a
finite number of discrete elements. And, these elements
were connected using multiple normal and shear springs.

Figure 1 Introduction of proposed method to theengineered (left) and non-engineered infilled frames.
(Right) (Colour code:
Grey = Reinforced Concrete, Pink = Soft Material, White= Weak Mortar and Red = Masonry Wall)
Table 1 Reduction in lateral drift of infill wall at 10mm
lateral load
Frame Type%drift reduction in
infill wall
Engineered Frame
84%
Non
Engineered 86%
Frame
3. CONCLUSION
The introduction of interface material reduces thetransfer of lateral forces on infill wall. The above Table 1
and Figure 1 demonstrate a significant reduction in lateraldisplacement and damage in infill wall, respectively. This
reduction prevents in-plane damages and renders
deformability to RC frame. The mitigation of in-plane
damages will intact infill wall strength and allow it tocarry the vertical load during failure of the reinforced
concrete structural members.
REFERENCES[1] R. Villaverde, “Seismic design of secondary structures:
State of the art,” J. Struct. Eng., vol. 123, no. 8, pp. 1011–
1019, 1997.
[2] H. Tagel-Din, “A new efficient method for nonlinear,large deformation and collapse analysis of structures,” Ph.
D. thesis, Civ. Eng. Dept., Univ. Tokyo, 1998.

― 70 ―

INFLUENCES OF ENTRAPPED AIR VOIDS ON THE MECHANICAL AND STRENGTHENING
CHARACTERISTICS OF SPRAYED FIBER-REINFORCED POLYMERIC (SFRP) COMPOSITES:MICROMECHANICS-BASED PARAMETRIC ANALYSIS
B.J. YANG1, HAMMAD, R. KHALID1, and H.K. LEE1
Department of Civil and Environmental Engineering,
Korea Advanced Institute of Science and Technology (KAIST),
291 Daehak-ro, Yuseong-gu, Daejeon 305-338, South Korea
Correspond to Prof. H.K. LEE (leeh@kaist.ac.kr)
1

Keywords: Micromechanics-based damage model, air voids, strengthening and mechanical characteristics, SFRP
composites
1. INTRODUCTION
In this study, the structural behavior of reinforced
concrete strengthened by the sprayed fiber-reinforced
polymeric composites (SFRP) composites considering
the influences of entrapped air voids is computationallyinvestigated.
A micromechanical constitutive damage model, which
accounts for the entrapped air void [1,2], is implemented
into finite element (FE) code ABAQUS [3], and theoverall strengthening performance of the system ispredicted based on the method proposed by [4].
2. NUMERICAL RESULTS
A numerical study considering the initial air voids is
performed to estimate the mechanical properties of SFRP
composites [4,5]. It shows that increase in the volume
fraction of initial air voids weakens the overall properties
of SFRP composites. In addition, the simulation result
reveals that higher volume fraction and aspect ratio offiber lead to stiffer stress-strain responses of thecomposite.
A comparison between the present numericalsimulations and experimental results from four-point
loading tests is also discussed in this paper [4,5]. The
predicted results are shown to be similar to theexperimental results, showing applicability of the model
to simulate the behavior of concrete beams strengthened
by SFRP composites (See, Figure 1). Details on theinfluences of entrapped air voids on the mechanical andstrengthening characteristics of SFRP composites will be
discussed.

3. CONCLUSIONSIn this study, a parametric analysis based on
micromechanics and FEM is conducted to evaluate theinfluences of entrapped air voids on the effective
behavior of SFRP composites. In the future, the proposed
approach would be further compared with additional
experimental data for the validation purpose.
ACKNOWLEDGEMENT
This research was supported by a grant (13SCIPA01)
from Smart Civil Infrastructure Research Program
funded by Ministry of Land, Infrastructure and Transport(MOLIT) of Korea government and Korea Agency for
Infrastructure Technology Advancement (KAIA).
REFERENCES[1] A. A., Karabutov and N. B. Podymova (2013).
“Quantitative analysis of the influence of voids anddelaminations on acoustic attenuation in CFRP
composites by the laser-ultrasonic spectroscopy method,”
Composites Part B: Engineering, Vol. 47, pp. 230-237.
[2] L. C. Hollaway (1994). “Handbook of polymer
composites for engineers,” The British Plastics
Federation.
[3] ABAQUS Version 6.5. Hibbitt, Karlsson & Sorensen
Inc, 2005.
[4] B. J. Yang, S. K. Ha, S. H. Pyo and H. K. Lee (2014).
“Mechanicalcharacteristics
andstrengthening
effectiveness of random-chopped FRP composites
containing air voids,” Composites Part B: Engineering,Vol. 62, pp. 159-166.
[5] Yang, B.J., Ha, S.K., and Lee, H.K. (2012).
“Micromechanics-based analysis on tensile behavior of
the sprayed FRP composites with chopped glass fibers,”
Journal of the Computational Structural Engineering
Institute of Korea, Vol. 25, No. 3, pp. 211-217. (inKorean).

Figure 1. The von-Mises stress of a reinforced concrete(RC) beam strengthen with SFRP composites― 71 ―

EQUATION FOR PREDICTIVE HORIZONTAL RESPONSE SPECTARA IN KOREA
D.H. SHIN1, S.J. HONG2, and H.-J. KIM3
Ph.D Student, Department of Architectural Engineering, University of Seoul, Seoul, Korea,
Graduate Student, Department of Architectural Engineering, University of Seoul, Seoul, Korea,
3
Associate Professor, Department of Architectural Engineering, University of Seoul, Seoul, Korea,
Correspond to Prof. H.-J. KIM (hyungjoonkim@uos.ac.kr)
1

2

Keywords: response spectra, nonlinear optimization, regression analysis, predictive equation
1. INTRODUCTION
Considering the characteristics of earthquake records
measured at Korea, this study empirically proposes a
predictive equation for horizontal response spectra thatcan be applicable to structural seismic design.
2. DERIVIATION OF THE ESTIMATINGEQUATION FOR RESPONSE SPECTRA
Total 1,404 ground motion records from 20
earthquakes in Korea are organized from 120 nationwide
earthquake stations of Korea institute of geoscience andmineral resources (KIGAM) and Korea meteorologicaladministration (KMA). Using this database, thecharacteristics of spectral shapes are investigated anddescribe in terms of shape parameters such as peakintensity (I), predominant period (Tsp), and width ofspectral shape (S) with respect to three independent
parameters such as magnitude (M), distance (R), and theaverage shear wave velocity in upper 30m (VS30).
As shown in Figure 1, the shapes of response spectra
are formulized as an approximate spectral shape function(SAnorm(T)) which is empirically determined from the
combination of the function F1(T) and F2(T) (Graizer andKalkan, 2009). Using the function, SAnorm(T), theapproximate spectral shapes of individual groundmotions can be accurately estimated by adjustment
procedure of shape parameters. In this study, theparameters representing the effects of basin depth, Dspand  are 1.75 and 1.25 from best fit.
This study aims to fit into the 3D response spectra
plotted in increasing order of distances using a single
continuous equation which can be used for all possible
earthquakes. To do this, the correlations of dependent
shape parameters are expressed as the sub-functions withrespect to three independent variables (M, R, VS30). In
order to find the best representative sub-function forshape parameters, this study performs the nonlinear
optimization coupled with regression analyses. This
procedure is repeatedly carried out for entire earthquake
data sets. From this result, estimator coefficients (m1~m4,
t1~t4 and s1~s3) minimizing difference between prediction
results and measured data are obtained. Finally, subfunctions of shapes parameters are substituted to baseline
equation, and response spectra are then predicted byanchoring the approximate spectral shape function topeak ground acceleration (PGA). Figure 2 comparatively
presents measured and predicted average response
spectra of several earthquake records.

Figure 1. Predictive equation for estimating response
spectra in Korea

Figure 2. Application of the predictive equationcomparison with major earthquake events in Korea
3. CONCLUSIONSThis paper makes an effort to propose empiricalequation for predicting response spectra in Korea. Basedon the proposed equation, predicted response spectra ofmajor earthquake events are well-matched withmeasured response spectra. Especially, peak spectral
acceleration, predominant period and slope of post-peak
decay are approximately corresponded to those ofmeasure response spectra.
ACKNOWLEDGEMENT
This research was supported by a grant from the
Research for improvement plan of existing public
structures in Korea funded by the Ministry of Public
Safety and Security.
REFERENCES[1] Graizer, V and Kalkan, E., Prediction of Spectral
Acceleration Response Ordinates Based on PGA
Attenuation, Earthquake Spectra, 2009, 25(1), pp. 39-69.

― 72 ―FINITE ELEMENT ANALYSIS OF REINFORCED CONCRETE BEAMS REPAIRED WITHULTRA-HIGH PERFORMANCE FIBER REINFORCED CONCRETE (UHPFRC)
M. SAFDAR1, T. MATSUMAOTO2, and K.KAKUMA3
Student, Graduate school of Engineering, Hokkaido University, Sapporo, Japan,
Associate Professor, Faculty of Engineering, Hokkaido University, Sapporo, Japan,
3
Researcher, Civil Engineering Research Institute for Cold Region, Sapporo, Japan,
Correspond to M. Safdar (safder_ali90@yahoo.com)
1

2

Keywords: ultra-high performance fiber reinforced concrete, finite element analysis, flexural strength, strengthening

2. NUMERICAL ANALYSIS
In finite element analysis, ordinary concrete is modeled
as an isotropic material capable of cracking in tension andcrushing in compression using Buyukozturk concretemodel. UHPFRC is also modeled as an isotropic material
with parabolic constitutive relation for compression andbilinear constitutive relation in tension. Truss elements areused for modeling of steel reinforcement (rebars andstirrups) with the common node points defined for eachrebar element and concrete solids. A bilinear constitutive
relation and von Mises failure criteria are employed forsteel modeling.
3. RESULTS
Figure 1 shows the load-deflection curves of repaired
beams and a reference beam (B-0) obtained by experimentand finite element analysis. Experimental and analyticalresults indicate that the ultimate flexural strength of RC
beams repaired with UHPFRC in flexural tension andcompression zone increased, with the increase of
UHPFRC thickness. This is attributed to the high tensile
strength and strain hardening of UPHFRC. Thicker
UHPFRC layer leads to increase in stiffness and cracking
load, thus, delaying the formation of localized macrocracks results in improving the protection function underservice conditions.

160

Compression zone repair

Load / KN

12080
EXP-B-0
EXP-BU-20
EXP-BU-40
EXP-BU-60

40

0

0

40

240

FEM-B-0
FEM-BU-20
FEM-BU-40
FEM-BU-60

80
Deflection / mm

120

160

Tension zone repair

200
160

Load / KN

1. INTRODUCTION
Over the past few decades, premature deterioration ofreinforced concrete structures exposed to severe
environmental conditions and mechanical loading has
become a serious problem. Previous studies have shownthat the use of ultra-high performance fiber reinforced
concrete (UHPFRC) improves the structural response andextends the durability of concrete structures.
In this study, the flexural behavior of reinforced
concrete beams retrofitted with UHPFRC is investigated,and experimental results are compared with 3-D finite
element analysis. The experiments were performed on
reinforced concrete beams repaired in tension andcompression zone, with UHPFRC of varying thicknesses.
The flexural strength of repaired beams was investigatedusing four-point bending test and compared with that ofreference beam without repair. After that, three
dimensional analysis was carried out by using a computer
software i.e. MSC/Marc to investigate the flexural
capacity and failure modes of all beams.

120
EXP-B-0
EXP-BL-20
EXP-BL-40
EXP-BL-60

80
40
0

0

40

80
Deflection / mm

FEM-B-0
FEM-BL-20
FEM-BL-40
FEM-BL-60
120

160

Figure 1. Comparison between the repaired beams andthe reference beam
4. CONCLUSIONSIn this study, the use of UHPFRC as restorative material
cross-section leads to higher stiffness and increased
bending capacity compared to previously repaired
thickness. In future based on the results of this study andadditional experiments, the structure performance of
UHPFRC as a strengthening material can be utilized.
REFERENCES[1] A.J. Wolanski, Flexural Behavior of Reinforced and
Prestressed Concrete Beams Using Finite Element
Analysis, Master’s Thesis, Marquette University,
Milwaukee, Wisconsin, May 2004.
[2] E. Brühwiler, E. Denarie, Rehabilitation andStrengthening of Concrete Structures using Ultra-High
Performance Fiber Reinforced Concrete, Concrete Repair,
Rehabilitation and Retrofitting III, 2012, pp.72-79.

― 73 ―RESPONSE ESTIMATION METHOD OF REINFORCED CONCRETE BUILDINGSDUE TO WATERBORNE DEBRIS IMPACT LOADS
T. ASAI1, K. MATSUKAWA2, H. CHOI2, and Y. NAKANO3
Graduate Student, Department of Architecture, The University of Tokyo, Tokyo, Japan,
Research Associate, Institute of Industrial Science, The University of Tokyo, Tokyo, Japan,
3
Professor, Institute of Industrial Science, The University of Tokyo, Tokyo, Japan,
Correspond to T. ASAI (asai-t@iis.u-tokyo.ac.jp)
1

2

Keywords: Impact Load, Tsunami, Tsunami Evacuation Building, Elastic Response Estimation
After the 2011 Great East Japan Earthquake, tsunami
loads were evaluated based on the damage observations [1],and a structural design code for tsunami evacuationbuildings was established [2]. However, impact loads due
to waterborne debris were not considered quantitatively
in the code. Waterborne debris were reported to cause
damage to buildings [3] (Figs. 1 and 2), and design
methods against the debris impact loads are therefore
currently in urgent need for designing safe tsunami
evacuation buildings.

The defined loads are secondly imposed to a particular
six story reinforced concrete building designed to an
expected tsunami wave load, and its elastic responses arecomputed by applying modal analyses. Because thebuilding response due to an impact load is significantly
affected by the impact duration, the response
characteristics are parametrically analyzed to different
impact durations. It is then found that the time, either
before or after the termination of the loading, when thestory drift reaches its maximum value depends on theimpact duration. Two simplified approaches are therefore
proposed considering impact duration to roughly
estimate the maximum story drifts due to the impact
loads.
As a result, the estimated maximum story drifts arefound in a good agreement with those computed bymodal analyses (Fig. 3).

Story drift ratio (%)

0.15

Figure 1. Building collapsed by ship collision(Kamaishi city, 2011) [3]Loaded floor : 6FL
Response floor: 2FL

0.10

0.05

0.00
0

0.1 0.2 0.3 0.4 0.5
Impact duration (s)Computed by modal analyses

Estimated

Figure 3. Example of estimated maximum story drifts

Figure 2. Stranded cargo ship collided with building(Kamaishi city, 2011)
In the previous researches, impact loads due towaterborne debris have been examined throughexperiments [4], but the building responses were notfocused in the discussions. To understand the building
responses due to the debris impact loads, in this paper,
impact loads due to ships and shipping containers are
firstly defined as a rectangular pulse, based on the studies
to evaluate collision strength of a ship’s bow [5] and thoseto obtain a time history of an impact load due to ashipping container [6].

REFERENCESG[1] T. Asai, T. Tateno, Y. Nakano, H. Fukuyama, K.
Fujima, Y. Haga, T. Sugano, T. Okada, Journal of
Structural Engineering Vol. 58B (2012) 97-104.
[2] Ministry of Land, Infrastructure, Transport andTourism, Notification 1318, Tokyo, 2011[3] Port and Airport Research Institute, Japan, Technical
Note of the Port and Airport Research Institute, No.1231,
Yokosuka, 2011[4] N. Mizutani, A. Usami, T. Koike, Annual Journal of
Civil Engineering in the Ocean, No. 23 (2007) 63-68.
[5] H. Nagasawa, K. Arita, M. Tani, S. Oka, Journal of
the Japan Society of Naval Architects and Ocean
Engineers, 142 (1977) 323-332.
[6] P. P. Aghl, C. J. Naito, H. R. Riggs, Journal of
Structural Engineering, 140 (5) (2014) 04014006.

― 74 ―

CODE2015

Technical Program – Day 2

Session 8

FULL-WAVEFORM INVERSION FOR MATERIAL PROFILE RECONSTRUCTION IN
PML-TRUNCATED DOMAINS
J.W. KANG1, H.D. Lee2, and B. KIM2
Assistant Professor, Department of Civil Engineering, Hongik University, Seoul, Korea,
Graduate Research Assistant, Department of Civil Engineering, Hongik University, Seoul, Korea,
Correspond to Prof. J.W. KANG (jwkang@hongik.ac.kr)
1

2

Keywords: Full-waveform inversion, Shear wave velocity, Perfectly matched layers (PMLs), PDE-constrained
optimization, KKT conditions, Regularization, Source frequency continuation
1.
ELASTODYNAMIC
FULL-WAVEFORM
INVERSION
In this study, a full-waveform inversion method formaterial profile reconstruction in heterogeneous semiinfinite domains was developed. In particular, themethod aims at imaging the spatial variation of shear
wave velocities, directly in the time-domain, from scantsurficial measurements of the domain’s response toprescribed dynamic excitation.
To truncate the semi-infinite extent of the physicaldomains in numerical analysis, perfectly-matched-layer(PML) wave-absorbing boundaries were introduced.
Within this framework, a mixed displacement-stressformulation based on unsplit-field PML was developed
for transient elastic wave simulations in heterogeneoussemi-infinite solid domains. The developed PML method
can be implemented by a mixed finite element method,and is relatively straightforward when compared withsplit-field PML methods. It also bypasses the need forcomplicated convolutional time integration which arises
when recent displacement-based formulations are used[1].
To formulate the full-waveform inversion of elastic
waves, a PDE-constrained optimization approach wasused, which lead to a classic KKT (Karush-Kuhn-Tucker)
conditions comprising an initial-value state, a final-value
adjoint, and a time-invariant control problem. The shear
wave velocity profile of medium was iteratively updated
by solving the KKT conditions with a reduced spaceapproach. To narrow the feasibility space and mitigate
the inherent solution multiplicity of the inverse problem,
Tikhonov (TN) and Total Variation (TV) regularization
schemes were used with a regularization factor continuation
algorithm. A source frequency continuation scheme wasalso utilized to make successive iterates remain within
the basin of attraction of the global minimum.
2. NUMERICAL RESULTS
Consider a 60m×30m truncated heterogeneous halfplane, where PMLs are introduced on the sides andbottom of the truncated domain, as shown in Figure 1(a).
The domain has 3 horizontal layers of shear wave
velocities 100 m/s, 115 m/s, and 130 m/s from surface tobottom. The density is 2000 kg/m� for all three layers.
Figure 1(b) exhibits the target velocity profile.
Figure 2 shows the reconstructed shear wave velocity
profiles using both TN and TV regularization schemes.
As shown in the figure, TV regularization captured the

sharp interfaces of medium better than TN scheme.

(a)

(b)

Figure 1 (a) Configuration of a PML-truncated heterogeneous
domain; (b) Target shear wave velocity profile(a) TN regularization(b) TV regularization

Figure 2 Inverted shear wave velocity profiles
In Figure 3, the variation of misfit, which is thedifference between calculated and measured responses, isplotted versus the number of iterations. The misfit
decreases monotonically as the iteration progresses.

Figure 3 Misfit versus number of inversion iterations
3. CONCLUSIONS
A full-waveform inversion method for reconstructing
the material profile of heterogeneous semi-infinite elastic
media was developed based on a PDE-constrained
optimization approach. Using the developed method, alayered shear wave velocity profile could bereconstructed with the aid of regularization schemes.
REFERENCES[1] U. Basu, A.K. Chopra, Int. J. Numer. Methods Eng.
59 (2004) 1039–1074.

― 75 ―

STUDY ON PORTABLE CELLPHONE DETECTOR FOR EARTHQUAKE
S.Y. Tang1, X.M. Shu2, S.F. Shen3, J. Hu1
Ph.D Candidate, Department of Engineering Physics, Tsinghua University, Beijing, China,
Associate Professor, Department of Engineering Physics, Tsinghua University, Beijing, China,
3
Professor, Department of Engineering Physics, Tsinghua University, Beijing, China,
Correspond to Prof. X.M. Shu (shuxm@tsinghua.edu.cn) and Prof. S.F. Shen(shensf@tsinghua.edu.cn)
1

2

Keywords: Fake base station, Cellphone detection, earthquake, Disaster relief.
1. INTRODUCTION
Mobile phone is one of the best beacon for detecting
survivor in disaster situation. A survey form project “ILOV” shows that there are about 80% victims carry theirmobile phone with them [1]. For the rescue in earthquake,a cellphone detecting and positioning system is proposedin this paper. This system is designed for searching
survivors indirectly by detecting and positioning
cellphones.
2. SYSTEM ARCHITECTURE
The cellphone detector is based on software definedradio (SDR) technology and the open software OpenBTS.

Figure 1. OpenBTS system diagram [2]OpenBTS builds Global System for Mobile
Communication (GSM) cellphone network, allow
cellphone to register in it and get International Mobile
Subscriber Identification number (IMSI) for individual
recognition. By initiate a voice call procedure, OpenBTS
can get transmit power of a cellphone and Receive Signal
Strength (RSS) of the system. By working with different
types of antenna, the system can get cellphone’s position
with triangulation algorithm of RSS and Angle of Arrival(AOA).
Hardware of the system include an integrated backpack
and an antenna. Figure 2 is the whole hardware system. A
series of component include Universal Software Radio
Peripheral (USRP), power amplifier, low noise amplifierand duplexer are in the backpack. The laptop run
OpenBTS and positioning algorithm. An omnidirectional
antenna, Yagi-Uda antenna or sector antenna can beapplied in this system.
The max rate of transmit power is 20W. With sector
antenna plugged into backpack, the system can detect
cellphones in extremely wide rage in disaster situation andrun AOA positioning algorithm. Omnidirectional antenna

is used to detect cellphone in all direction and run RSS
positioning algorithm. Yagi-Uda antenna is used for

poximity sensing.
Figure 2. Hardware system
3. CONCLUSIONSWith several experiments, some key features of thissystem are proven. (1) The system allow all cellphonewhich support GSM protocol to register in, and accept
IMSI in the same time; (2) By start a voice call procedure
with a specific cellphone, the system can continuously
accept cellphone transmit power and RSS; (3) When thesystem is running with directional antenna, the angle ofsignal arrival can be reached by the correlation between
antenna direction and RSS; (4) Working with sector
antenna in full transmit power, cellphones in 1.2km away
are able to register in.
The results shown above, indicate that the system issuitable to work in huge disaster situation, especiallyeffective in no mobile service network case. The
cellphone detecting and positioning system can detect
mobile phones in more than 4km2 range in short time,acquire rough number of victims and evaluate the demandof rescue resource.
REFERENCES[1] L. Chen, T. Ostertag, M. Loschonsky, L.M. Reindl,2010 IEEE International Conference on WirelessCommunications, Networking and Information Security,
IEEE, New Jersey, 2010, pp. 252–256.
[2] K. Vasudeva, B.S. Ciftler, A. Altamar, I. Guvenc,
2014 IEEE 15th Annual Wireless and Microwave
Technology Conference, IEEE, New Jersey, 2014, pp. 1–
6― 76 ―

CRACK PROPAGATION SIMULATION BY THE EXTENDED CONSECUTIVEINTERPOLATION FINITE ELEMENT METHOD (XCQ4)
Z.Y. KANG1, T.Q. BUI2, and S. HIROSE3
Ph.D student, Department of Mechanical and Environmental Informatics, Tokyo Institute of Technology, Tokyo, Japan,
2
Researcher, Department of Mechanical and Environmental Informatics, Tokyo Institute of Technology, Tokyo, Japan,
3
Professor, Department of Mechanical and Environmental Informatics, Tokyo Institute of Technology, Tokyo, Japan,
Correspond to Ph.D student Z.Y. Kang (kang.z.aa@m.titech.ac.jp)

1

Keywords: Fracture, XFEM, consecutive-interpolation, crack growth, stress intensity factor
1. INSTRCUTIONSThe extended finite element method [1] is originallydeveloped for crack growth without remeshing in termsof the partition of unity method. Currently, Bui et al. [2]developed a new consecutive-interpolation 4-nodequadrilateral finite element (CQ4) initially applied tostress analysis of 2D solids. Based on the CQ4 scheme
the extended consecutive-interpolation finite element islater developed, named as XCQ4, for stationary cracks in
2D solids [3]. The numerical results indicate that the
XCQ4 offers a better accuracy and higher convergence
rate compared with those derived from the classical
XFEM using 4-node quadrilateral element (XQ4). The
underlying idea behind the proposed XCQ4 is that theapproximation functions constructed based on the CIP
contain both nodal values and averaged nodal gradientsas interpolation conditions.
In this work, the new element (XCQ4) is further
extended to predict the crack growth paths in 2D linear
elastic solids. We integrate the interaction integral interms of J-domain to estimate the fracture parameters.
2. FORMULATION OF THE XCQ4 ELEMENT
The enriched displacement approximation using the
new XCQ4 is executed in a similar way with XQ4 asu h ( x ) = ∑ N%i ( x )u i +
i∈I s

+

∑ N% j ( x )  H ( x ) − H ( x j ) a j

j∈J cut(1)

4∑ N%k ( x ) ∑  F α ( x ) − F α ( x k ) bαk
α =1k∈K tip

cut

where I , J , K tip denotes the total set, Heavisideenriched set and asymptotic enriched set respectively. In
addition the XCQ4 shape function is defined as Eq. (2).
nodex (x=i,j,k,m) denotes the surrounding nodes ofconsidered quadrilateral element.
s

N% f = φi N [fi ] + φix N [fi,]x + φiy N [fi,]y + φ j N [f j ] + φ jx N [f j, ]x + φ jy N [f j, ]y
1 4 4 44 2 4 4 4 43
1 4 4 4 4 2 4 4 4 43
node inode j+ φk N [fk ] + φkx N [fk, ]x + φky N [fk, ]y + φm N [fm ] + φmx N [fm, x] + φmy N [fm, y]
1 4 4 4 4 2 4 4 4 43
1 4 4 4 42 4 4 4 43
node k(2)

node mIn terms of the maximum hoop stress criterion thecrack should initiate in a direction defined as Eq. (3)when the circumferential stress is maximum.
2

K 
1 K
θc = 2 tan −1 ( I ±  I  + 8)
4 K II K II (3)

3. NUMERICAL RESULTS
A cracked beam with three holes is adopted to assess

the robustness of XCQ4. Two test cases having different
values of a and b is computed as case I a = 1.5 b = 5.0 ;case II a = 1.0 b = 4.0 . In both cases the numericalresults show a desirable agreement with experimental
solutions reported by Bittencourt et al. [4].

(a)

(b)
(c)
Figure 1. (a) Configuration of cracked beam withthree holes; Comparison of predicted crack paths withexperimental solution (b) case I; (c) case II
4. CONCLUSIONS
A novel approach for modelling crack propagation was
proposed. The XCQ4 method is successfully validated as
the predicted crack growth paths agree well with reliable
experimental solutions.
REFERENCES[1] N. Moes, J. Dolbow, T. Belytschko, A finite element
method for crack growth without remeshing. Int. J.
Numer. Meth. Engng., 46 (1999) 131–150.
[2] T.Q. Bui, Q.V. Dam, Ch.zhang, D.N. Du, A
consecutive-interpolation quadrilateral element (CQ4):
Formulation and applications. Finite Elem. Anal. Des.,84 (2014) 14-31.
[3] Z.Y. Kang, T.Q. Bui, D. Nguyen, T. Saitoh, S. Hirose,
An extended consecutive-interpolation quadrilateralelement (XCQ4) applied to linear elastic fracture
mechanics. Acta Mech.
[4] T.N. Bittencourt, P.A. Wawrzynek, A.R. Ingraffea,
Quasi-automatic simulation of crack propagation for 2d
LEFM problems. Engineering Fracture Mechanics., 55(1996) 321-334.

― 77 ―

SIMULATION OF FAILURE PROCESS OF REINFORCED CONCRETE BEAM COLUMNJOINT WITH MECHANICAL ANCHORAGES BY 3D RBSM
L.EDDY1 and K. NAGAI2
Graduate Student, Department of Civil Engineering, The University of Tokyo, Tokyo, Japan,
Associate Professor, International Center for Urban Safety Engineering, Institute of Industrial Science, The Universityof Tokyo, Tokyo, Japan
Correspond to L.Eddy (eddy@iis.u-tokyo.ac.jp)
1

2

Keywords: Mechanical Anchorages, 3D RBSM, Knee Beam Column Joint, Failure Process

Load (kN)

Simulation
Def.×3

100

Experiment50
Simulation
Experiment

0

0
10
20
30
40
Displacement (× 10-3 rad)Figure 2. Load-displacement
Figurerelationships
150 Bond works Diagonal cracks
1.

2.

100

Load (kN)

1. BACKGROUND
The use of mechanical anchorages is still limited when
they are placed near the surface of the beam column joint
since the internal stress and cracking pattern have notbeen well understood. A meso-scale analysis by 3D
discrete method, called 3D RBSM, is conducted, because
the 3D shape of a reinforcement bar is modeled directly
including ribs of a reinforcement bar, cracks occur as theresults of discontinuous of concrete and interaction
between concrete and a reinforcement bar at meso-scalelevel, and cracks can be simulated directly. In this study,
by investigating the internal stress and cracking patternof simulation results, the failure process of the beamcolumn joint with mechanical anchorages is revealed
through the comparison with the experimental results.
In RBSM, proposed by Kawai et al. (1978), a 3D
reinforced concrete member is meshed into rigid bodies.
The simulation system, developed by Nagai et al. (2005),
is used.

Diagonal cracks
open

Cracks propagate
to the surface

50

3.Surfacecracks4.

3.
2. NUMERICAL MODEL
Simulation is conducted for the experiment of beamcolumn joint with mechanical anchorages done by
Yoshimura et al. (2012). 756642 elements are used tomodel the beam column joint (Figure 1). The same
dimension and material properties as experiment aremodeled.
3. SIMULATION RESULTS
Figure 2 shows the load displacement relationship ofsimulation result which is roughly the same as that ofexperimental result with 10% difference of capacity, i.e.

430

300

Column

1915

D10

Rigid
plate
480 480

Pin
Rigid plate

D19
D22
D10

Loaded Joint Detail
Pin

Figure 1. Detail of numerical model (Units: mm)0010Drift Angle (× 10-3)

20

Figure 4. Failure process based on the simulation result
85 kN and 94.5 kN in case of simulation and experiment,
respectively. Figure 3 shows the surface cracks of thenumerical model which are also roughly the same as
those of experimental specimen. Based on the study ofthe internal stress and cracking pattern of simulation
result, the failure process is revealed, shown on Figure 4.
4. CONCLUSIONS
The failure process of the beam column joint withmechanical anchorages has been revealed through thestudy of internal stress and cracking pattern of 3D RBSM.
REFERENCES[1] T. Kawai., New Discrete Models and Their
Application to seismic Response Analysis of Structure,
Nuclar Engineering and Design, 1978, 48, pp. 207-229.
[2] K. Nagai, Y. Sato, and T. Ueda, Mesoscopic
Simulation of Failure of Mortar and Concrete by 3D
RBSM, J. Adv. Conc. Technol., 2005, 3(3), pp.385-402.
[3] M.Yoshimura, T. Kiyohara, A. Tasai, and K.
Kusunoki, Experimental Study of Performance
Improvement of Knee Joints with Mechanical Anchorage,
Proceedings of JCI, 2012, 34(2). (in Japanese)― 78 ―FINITE ELEMENT ANALYSIS OF BURIED PIPELINE ON FORST REGION CONSIDERINGFROST HEAVE
K.W. KIM1, M.K. KIM2, W.Y. CHO3, and K.Y. YEAU3
Ph.D candidate, Department of Civil and Environmental Engineering, Yonsei University, Seoul, Republic of Korea,
2
Professor, School of Civil and Environmental Engineering, Yonsei University, Seoul, Republic of Korea,
3
Senior Researcher, Energy Infrastructure Research Group, Steel Solution Center, POSCO, Incheon, Republic of Korea,
Correspond to Ph.D Candidate. K.W. KIM (kyu_won@yonsei.ac.kr)
1

Keywords: frost heave, soil freezing, porosity rate, buried pipeline
1. INTRODUCTION
Due to the growing demand for oil and gas in theworld, there is an increasing interest to explore new areasfor petroleum production. The Arctic Circle has 13% of
the undiscovered oil in the world. Oil and gas aretransported over long distances in pipelines, and thepipeline crosses regions of permafrost or seasonal frost.
Ground-induced actions due to frost heaving and thawing
in regions of seasonal frost are responsible for the
majority of damages in oil and gas buried steel pipelines.
Therefore, researches of pipelines crossing frozen soils
have been performed and several models describing theheaving process have been developed in the past.
The porosity rate model is indicative of ice growth. Ice
growth is dependent on the temperature gradient and thestress state in the freezing soil. The porosity rate model
was introduced by Michalowski and Zhu [1], and themodel was implemented in a 2-D Problems. When soils
are frozen by air or inner pipeline temperature, ices madethe adhesive force between the pipeline and the soil. The
3-D numerical model is necessary to consider thisphenomenon. Therefore, this study performs numericalanalysis for the 3-D pipeline using the porosity rate
model. We compared the numerical results between the
1-D model [1] and the 3-D model.
2. NUMERICAL ANALYSISThe numerical computations are performed using thecommercial program ABAQUS. Some assumptions areapplied in this study. The heat transfer takes place byconduction only. The soil is fully saturated. The thermal
conductivity of soil is isotropic. The behavior of the soil
to loads is elastic, but the elastic properties depend on thetemperature. The total strain increment consists of boththe elastic strain increment and the growth increment due
to the change in porosity. The thermal constitutive
behavior of the soil and the porosity rate model aredefined and integrated into ABAQUS using the user
subroutines UMATHT and UEXPAN, respectively.
Two types of problems are performed in this study: 1D freezing (1-D heat flow) for validation of the model,
and 3-D pipeline analysis. The analysis is a fully coupled
thermal-displacement analysis.
The pipeline is buried 1.72m below the top surface(Figure 1). Pipeline is composed of steel. The diameterof the pipe is 0.9m, and the thickness of the pipe is
8.5mm. The pipeline is not insulated. The width of themodel is 12.0m and the height is 8.0m. The size of theanalysis model was formed large enough to determine
the temperature distribution due to the change of air andinner pipeline temperature. Soil skeleton, water, ice
thermal properties that used in this study are represented
in Table 1. The total simulation time was 60days.

Figure 1. Finite element modelling
Table 1. Thermal properties of soil constituents
Mass heat
ThermalDensitycapacity
conductivity(kg/m3)(J/kg℃)(W/m3℃)
Soil
2620
900
2.92
skeleton
Water
100
4180
0.56
Ice

917

2100

2.24

3. CONCLUSIONSThis study was extended to 3-D model from 2-D
porosity rate model to consider the effect between thepipeline and ground. As a result of numerical analysis,
the presented model in this study confirmed thatindicated the 3-D behavior of the frozen ground. But
many experiments and data of field are still needed toindicate the realistic behavior of the frozen ground. The
model of this study is expected that be able to consider
the behavior of discontinuous frozen ground or the pipesoil interaction.
REFERENCES[1] R.L. Michalowski, M. Zhu, Frost Heave Modelling
Using Porosity Rate Function, Int. J. Numer. Anal. Mech.
Geomech., 30 (2006) 703-722.

― 79 ―

STRAUTO : PARAMETRIC STRUCTURAL MODELLING AND OPTIMIZATION SYSTEM
C.K. KIM1, S.E. LEE2
Professor, Dept. of Architectural Eng., Dankook University, Yongin-si, Korea,
Professor, School of Information & Communication, Baekseok University, Cheonan-si, Korea
Correspond to Prof. C.K.KIM (cheekim@dankook.ac.kr)
1

2

Keywords: StrAuto, Parametric Design, Computational Optimization, Distributed Parallel Solver System
1. INTRODUCTION OF STRAUTO
StrAuto is a parametric structural optimization tool
developed by our team, a plug-in for Grasshopper runningon Rhinoceros. It uses SAP2000 and ETABS APIs for
operating real-time iterative analysis, controlling theexternal solvers and bringing back the results in real-timefor evaluation. Thus, selecting an optimum design isenabled by utilizing StrAuto’s iterative analysis function,
which checks a number of automatically generated
alternatives and brings back their analysis results in realtime for evaluation.
A key feature of StrAuto’s design optimization function
lies in the iterative data exchange between auto-generated
geometric models and external analysis solvers. It is theAPIs that interconnect parametric analysis models,geometries of which can change in real time, with external
analysis solvers. StrAuto interacts with the solvers tochange the model structure’s geometry or properties,
manages analysis, and evaluates alternatives byprocessing analysis results. Until the evaluated results
reach a specific goal, StrAuto creates new combinationsof parameters to develop the analysis model. By suchfeatures, thousands of design alternatives, which can takeover a week for modelling a single alternative if using
existing engineering platforms, can be created, analyzed,
and compared in just a few hours for locating optimal
solution.
In order to achieve the lowest cost possible, multiple
simulation runs should be performed and project-specificprocess should be built for finding the fittest structural
system or member section. Material strengths and thesizes of vertical structural members (columns and walls)
are usually examined for tall building cases.

Figure 1. 4 StrAuto’s Design Optimization Process
2. APPLICATIONSThe Blade reaches 300 meters with the total gross floorarea of 101,700 square meters. It is diamond-shaped inplan with its sharpest edges at the north-east and southwest corners of the site.
We examined alternative structural systems for the
tower, mega brace and tube wall system to replace theoriginal structure. We tested multiple steel sections withmultiple vertical zonings for each of the three structural
systems to determine an optimal system for the tower. Not
only did we provide the comparative data for material
quantities for each system, but we also provided optimized
design results for each of them. The tube wall system wasproved to satisfy the minimum cost.

Figure 2. Comparative Optimization for The Blade
3. CONCLUSIONSStrAuto’s functionality is clear and powerful and can beeffectively adapted to custom projects. StrAuto’s
capability to check numerous alternative cases based on
the design development plans can be greatly utilized tosuggest guidelines for the subsequent construction designphase. We expect to expand StrAuto’s functionality tosupport multiple types of project cases based on our trackrecord of delivering diverse major projects.
REFERENCES[1] Kim, Y. et al. (2012) Selection of Optimal Structural
System for Complex-shaped Super-tall Building
Structural Systems Using Parametric Design Technique.
Journal of the Architectural Institute of Korea: Structure& Construction, 28 (1), 77-84[2] Kim, Y. et al. (2009) A Case Study on the Field
Application of Parametric Design Technique for
Optimization of Outrigger System High-rise Building
Structures. Journal of the Architectural Institute of Korea:
Structure & Construction, 28 (12), 45-52― 80 ―

CODE2015

Technical Program – Day 2

Session 9

NUMERICAL EVALUATION OF DAMAGE STATES AND SEISMIC CAPACITY OF
REINFORCED CONCRETE COLUMNS
C.K. GADAGAMMA1, and K. MEGURO2
Graduate Student, Department of Civil Engineering, The University of Tokyo(UTokyo), Tokyo, Japan,
Director/Professor, International Center for Urban Safety Engineering (ICUS), Institute of Industrial Science(IIS),
UTokyo , Tokyo, Japan,
1

2

Keywords: AEM, Reinforced Concrete Columns, PEER, Seismic Capacity
1. INTRODUCTION
To assess the seismic vulnerability of reinforced
concrete buildings, a damage database is being developed,
which requires an accurate and sophisticated numerical
tool. In this study, a discrete element approach based
numerical tool called Applied Element Method (AEM) [1]is used. To verify the capability of AEM to carry out
inelastic analysis and study the properties of plastic hinges
corresponding to different damage states, validation ofrectangular RC columns was required. Pacific Earthquake
Engineering Research (PEER) experimental database forstructural performance, which is a compilation of various
experiments conducted on rectangular reinforced columns
against monotonic cyclic loading was selected. This paper
evaluates visual damage states as provided in the databasewith respect to measurable parameters using numericalsimulation.
2. MODELLING OF RC COLUMNS USING AEM
In this study, test results were studied from theliterature [2], [3] and the most suitable ones have beenshortlisted for numerical validation. Figure 1 shows atypical AEM model of test specimen.
3. VALIDATION
Using AEM, the experimental results have been
numerically simulated and the corresponding results are
as shown in Figure 2. The necessity of further
investigation into the results for accuracy was desired.
This led to the examination of material models being usedin the existing AEM. The
uniaxial compression material
model being used in AEM isbased on Maekawa’s ElastoPlastic Fracture (EPF) model[4], which considers single post
peak behaviour of concrete incompression. Observing this,authors have considered thepossibility of size effects inconcrete during post peak
softening of concrete under
uniaxial compression. The
brittle nature of concrete as
Figure 1 AEM Modelseen in Figure 2, is mainly due
to size effects, which were of RC column by [1]initially
witnessed
andexperimentally verified the existence in 1990s [5]. Further,

the constitutive models for concrete were modified basedon parameters from fracture mechanics and validated andthe detailed analysis is being communicated as a part ofanother study. These combined models have been used to

Figure 2 Preliminary Results of [1] and [2]revaluate RC columns from PEER experimental database.
3. CONCLUSIONS
The results obtained from the simulations are inagreement with experimental results and with improved
accuracy due to the modification in the constitutive
material model of concrete subjected to uniaxial
compression. The visual damage states were evaluatedusing crack initiation and propagation. The changes inplastic hinge properties are evaluated by studying thespring properties for different damage states. Further, theseismic capacity of the columns has been discussed byassessing energy dissipation, stiffness and strength.
REFERENCES[1]H. Tagel-Din, “A new efficient method for nonlinear,large deformation and collapse analysis of structures.”東京大学, 1998.
[2]H. Tanaka, “Effect of lateral confining reinforcementon the ductile behaviour of reinforced concrete columns.”
1990.
[3]O. Bayrak and S. Sheikh, “Confinement steel
requirements for high strength concrete columns,” inProceedings of 11th World Conference on Earthquake
Engineering, 1996.
[4]K. Maekawa, H. Okamura, and A. Pimanmas, Nonlinear mechanics of reinforced concrete. CRC Press, 2003.
[5]J. B. M. Van Mier, “Strain-softening of Concrete under
Multiaxial Loading Conditions, PhD thesis,” pp. 1–244,
1984.

― 81 ―

DEVELOPING A NUMERICAL METHOD FOR FATIGUE LIFE PREDICTION
OF RC SLABS REINFORCED WITH PLAIN BARS
Ahmed Attia M. DRAR1 and Takashi MATSUMOTO2
1

Doctoral student, Graduate School of Engineering, Hokkaido University, Sapporo, Japan
2
Associate professor, Faculty of Engineering, Hokkaido University, Sapporo, Japan

Keywords: Fatigue, RC slabs, Bridging degradation, Moving load, Plain bars
1.

INTRODUCTION
Recently, there is an increasing request to estimate thefatigue behaviors of RC bridge slabs under a moving load.
Most of the previous studies focused on the modeling of
fatigue behaviors of RC slabs reinforced with deformed
bars. However, many RC slabs in use today are reinforced
with plain bars, and they are suffering from fatigue
damages. This study presents a numerical method tosimulate the fatigue behaviors of RC slabs reinforced withplain bars under a moving load.

3.

RESULTS AND DISCUSSION
Applying moving load leads to the cracked elements
propagation. By increasing the number of cycles, thebridging stress degradation of these cracked elements leadto a decrease in slab stiffness and an increase in slab
deformations. The center displacement evolution versus
the number of cycles for all slabs is compared with theexperimental results in Figure 2. The slab under a higher
moving load level shows a larger slab center displacement
than that under a lower moving load level due to the initialdisplacement at first cycle. In addition, at a higher moving
load level, the center displacement evolution shows alarger slope than that at a lower moving load level. The
reason is that the slab under a higher moving load level
deteriorated
more
quickly
than
Table 1. RC slabs details [4]SlabID

Movingload (kN)

Rebar type

P110
P150

110
150

Plain bar
SR235

Concretestrength
(MPa)
43.0
41.7

36.6

2650 mm3300 mmMoving load zone 300 × 2000 mm
D16 @ 260 mm

D16 @ 130 mm

Figure 1. RC slab geometry
15
Slab center displacement (mm)METHOD
An RC slab model of smeared crack elements is solved
using finite element method. A moving load wassimulated by moving a constant load level along thelongitudinal direction as shown in Figure 1. The cracked
elements are modified according to the bridging stressdegradation concept [1, 2]. This concept can be defined as
the reduction of transferred stress between crack surfaces
through aggregates under repetitive loading. The bondslip effect between a plain reinforcing bar and its
surrounding concrete is taken into consideration by adding
equivalent bond strain to plain bar strain [3]. The
numerical model is verified using the experimental databy Shakushiro et al. [4]. The details of the tested slabs areshown in Table 1.

190

D13 @ 230 mm

2.

P190

Exp.[4] Ana.
12
9

P110
P150
P190

6
3
0 0
10

1

10

2

3

4

5

10
10
10
10
Number of cycles (Cycles)

6

10

Figure 2. Center displacement evolutions
that at lower moving load level. The RC slab subjected toa higher moving load shows a shorter fatigue life than thatat a lower moving load. The experimental and numericalresults reveal similar values for fatigue life and center
displacement, indicating an acceptable agreement
between them.
4.

CONCLUSIONS
The fatigue analysis of RC slabs reinforced with plain
bars under moving load was conducted using thisnumerical method, and it shows a good agreement withthe experimental results. Bond-slip effect of plain bars
shall be taken into consideration for accurate fatigue life
prediction.
REFERENCES[1] V. C. Li, T. Matsumoto, Cement and concrete composites.
20(5) (1998) 339-351.
[2] Z. Jun, H. Stang, ACI Materials Journal. 95(1) (1998) 58-67.
[3] M. Dehestani, S. S. Mousavi, Construction and Building
Materials. 81 (2015) 284-290.
[4] K. Shakushiro, H. Mitamura, T. Watanabe, N. Kishi, Journalof Structural Engineering. 57(A) (2011) 1297-1304 (in Japanese).

― 82 ―FATIGUE ANALYSIS OF FROST DAMAGED RC SLABS SUBJECTED TO MOVING LOAD
BASED ON BRIDGING STRESS DEGRADATION CONCEPT
S. EDALO1 and T. MATSUMOTO2
Master student, Graduate School of Engineering, Hokkaido University, Sapporo, Japan,
2
Associate professor, Faculty of Engineering, Hokkaido University, Sapporo, Japan,
Correspond to S. EDALO (edaloserge@gmail.com)

1

Keywords: Fatigue, Frost Damage, RCSlab, Bridging Stress Degradation

2. METHOD
Finite element method is used to solve a slab model ofsmeared crack elements. The cracked elements are
modified according to the bridging stress degradation
concept by Li and Matsumoto [1]. Frost damaged concrete
stress-strain model by Hasan et al [2] is used to take into
account the frost damage. Giuffré-Menegotto-Pinto model[3] is employed to represent the stress-strain behaviour ofa rebar under repetitive load.
3. RC SLAB MODELINGThe studied slab is of dimension 2300 ×3000 × 180 mm
and is supported by steel I-beams along its width. Frostdamage is simulated in experimental study by reducing
the slab thickness by 10 mm rather than freezing-thawing
the slab [4].
On the other hand, in the present proposed finite
element method, the slab is divided into 5 layers, and the
top layer of 10 mm thickness is frost damaged. The
concrete and reinforcement properties of the slab are thesame as those used in the experimental study. A moving
fatigue load is applied along the slab length centerline asindicated in table 1.
Table 1. Applied load
Number of cycle N
1 ≤ N ≤ 100000
100000 < N ≤ 200000
200000 < N ≤ 300000
300000 < N ≤ 400000
400000 < N

Load (kN)
130
140
170
200
230

20
16

Deflection (mm)1. INTRODUCTION
Freeze-thaw cycling is one of the main causes ofreinforced concrete (RC) deck slab deterioration in cold
regions. This deterioration stems from water freezing in
concrete upon cooling, and its thawing upon subsequentheating. This process is accompanied by dimensionalchange and internal stress change, and leads to a
progressive accumulated physical damage in concrete.
Freeze-thaw cycling can result in severe damage.
Therefore, it is necessary to quantitatively evaluate therelationship between the level of frost damage and theremaining lifespan of RC deck slabs. This study aims topropose an experimentally verified numerical model for
predicting the fatigue life of frost damaged RC deck slabs.

EXP
EXP

ANA
ANA

12
8

With Frost Damage

4
0 0
10

Without Frost Damage

10

1

10

2

10

3

10

4

10

5

10

Number of cycles (Cycles)
Figure 1. Midspan displacement evolutions
4. RESULT
Figure 1 shows the midspan displacement evolutions
obtained by experiment and numerical analysis for theslab without frost damage and that with 10 mm frost
damage.
It comes out that the midspan displacement evolution
for the RC slab with frost damage shows a larger slope
than that without frost damage. A possible explanation forthis is the reduction of the slab flexural stiffness due to theconcrete stiffness deterioration caused by the frost damage.
The fatigue life of the RC slab with frost damage isshorter than that without frost damage consideration. Thisis due to the degradation of the concrete initial strengthand stiffness by the frost damage.
Numerical results show a good agreement withexperimental ones.
5. CONCLUSION
The good correlation between the experimental andnumerical result confirms the applicability of the proposedmodel to predict fatigue behaviour of frost damaged slab
under moving load.
REFERENCES[1] V. C. Li, T. Matsumoto, Cement and Concrete
Composites. 20(5) (1998) 339-351.
[2] M. Hasan, H. Okuyama, Y. Sato, T. Ueda, Journal ofAdvanced Concrete Technology. 2(1) (2004) 89-99.
[3] M. Menegotto, P.E. Pinto, IABSE reports of theworking commissions. 13(1973) 15-22.
[4] H. Mitamura, T. Satou, K. Honda, S. Matsui, Journalof Structural Engineering, JSCE. 55(A) (2009) 1420-1431(in Japanese).

― 83 ―

6

EFFECT OF THE BASE ISOLATION SYSTEM ON THE LNG STORAGE TANK USINGPROBABILISTIC METHOD
J.-S. KIM1, T.-H. LEE2, and T.-S. HAN3
Graduate student, Department of Civil and Environmental Engineering, Yonsei University, Seoul, Korea,
2
Associate Professor, Department of Civil Engineering, Konkuk University, Seoul, Korea,
3
Associate Professor, Department of Civil and Environmental Engineering, Yonsei University, Seoul, Korea,
Correspond to Prof. T.-S. HAN (tshan@yonsei.ac.kr)
1

Keywords: Base Isolation System, LNG Storage Tank, Probabilistic Method, First-order Second-moment Method,
OpenSees, Finite Element Method
1. INTRODUCTION
As the natural gas demand has been steadily increased,
the interest in liquefied natural gas storage tank is also
increased. When the LNG storage tank is designed, thelarge safety factor will be considered for earthquake load.
It is very hazardous that the liquefied natural gas isaccidentally spilled out or leaked by the earthquake.
Therefore, LNG storage tank will be protected fromrandom ground motions by seismic design. Among
various seismic design methods, using base isolation
system is one of the effective methods to reduce the lateralforce from earthquake. There are several uncertain
mechanical properties which affect the seismically
performance. To provide more accurate seismic analysis,
the probabilistic approach such as the first-order secondmoment (FOSM) method can be used. The objective of
this study is to investigate the effect of the base isolation
system on the LNG storage tank considering the
distributed properties of base isolation system.
2. PROBABILITY APPROACH FOR BASE
ISOLATION SYSTEM
The behavior of the base isolation system depends onthe mechanical and geometrical properties. Because each
property involves uncertainty and error, probability
approach would be beneficial to provide accurate seismic
analysis. The first-order second-moment (FOSM) method
is one of the most popular methods to determine thestochastic moments of a function with random input
variables [1]. In previous studies, random variables are
selected as structural properties such as concrete strength

Figure 1. Acceleration response spectra of the 20
different earthquakes. Bold line is average response
and initial stiffness or earthquake input profiles. The
property of base isolation device is one of the key factors

to seismic performance, and the FOSM method is applied
to analyze the seismic effect of the LNG storage tanks.
3.

MODELING PROCESSThe LNG storage tank (Lee et al., 2013) and the baseisolation system are modeled by OpenSees, a finite
element analysis software framework for earthquake
engineering simulations [2]. It is simply assumed that the
LNG storage tank consists of concrete outer tank, steel
inner tank, and liquid part as shown in Fig, 2. The ‘Tuningfork model’ is used to model the LNG storage tank. Thebase isolation system is located at the bottom of the tank,which is modeled using the ‘elastomericBearing’ element
provided by OpenSees. Using this model, nonlinear time
history analyses are conducted.
4.

CONCLUSIONS

Figure 2. Schematic of the LNG storage tank.
The effect of seismic isolation system is investigatedfrom nonlinear time history analyses when the propertiesof base isolator are considered as probabilistic variables.
Using the FOSM method, more accurate seismic analysis
can be provided and this result can be utilized to design
the optimal base isolation system.
REFERENCES[1] T.-H. Lee, K.M. Mosalam, Earthquake Engng. Struct.
Dyn. 34 (2005) 1719–1736.
[2] F. McKenna, G.L. Fenves, M.H. Scott, and B.
Jeremic, Open System for Earthquake Engineering
Simulation (OpenSees), Pacific Earthquake Engineering
Research Center, University of California, Berkeley, 2000.

― 84 ―MONITORING OF A FULL SCALE PRESTRESSD CONCRETE GIRDER
DURING CONSTRUCTION PROCESS
T.M. KIM1, D.H. KIM2, M.K. KIM3, and Y.M. LIM3
Ph.D. Candidate, Department of Civil & Environmental Engineering, Seoul, Korea,
Senior Researcher, Civil & Marine engineering research team , GS E&C, Seoul, Korea,
3
Professor, Department of Civil & Environmental Engineering, Seoul, Korea,
Correspond to T.M. KIM (taemin@yonsei.ac.kr)
1

2

Keywords: Monitoring, Full Scale, Prestressed Concrete (PSC), Construction Process
1. INTRODUCTION
Many other studies already have tried to measure thebehavior of girder using the fiber bragg grating (FBG)
sensors. However, it was hard to measure the behavior ofprestressed concrete (PSC) girder in the course of being
lifted and carried. During the construction process of
PSC girder, the structural system of it changes due to thechanges of loading and boundary conditions, so that themonitoring of these changes is important. Thus, thispaper deals with the method to monitor a full scale PSC
girder during construction process.
2. SENSOR INSTALLATION AND
MEASUREMENT
In order to monitor the deflection and strain, long
length FBG sensor was fabricated and applied to the PSC
girder. Full scale PSC girder of 60 m span was alsofabricated and three FBG sensors were attached at thebottom fiber of mid-span and both 1/4-spans. In the midspan, additional FBG senor was attached at the top fiber
to measure the neutral axis. In order ot install the longgauge FBG sensors, iron support devices were embedded
to the girder. Figure 1 shows the sensor attachment
locations and Figure 2 shows the actual installation of
FBG sensors.

process of lifting girder. In order to verify the effect ofchanged structural system in these process, occurredstrain and displacement were measured. The way toestimate of deflection was based on Euler–Bernoulli
beam theory.

Figure 3. Lifting Girder.
The girder was carried to laboratory and then four
point bending test was conducted. In this process, theway to estimation of deflection was verified bymeasuring of deflection using LVDT. Estimateddeflection from FBG sensors and measured deflection
from LVDT were compared.
3. CONCLUSIONSThis paper presents the way to estimate the strain of
the PSC girder using FBG sensors during theconstruction process. This method using long gauge FBG
sensor was able to measure the behaviors of girder
effectively.

Figure 1. Attachment Locations of FBG Sensors

Figure 2. Sensor installation.
After installation of sensors, the fabricated PSC girder
was lifted by crane and carried by trailer from the casting
bed to the structural laboratory. Figure 3 shows the actual

REFERENCES[1] W. Chung, S. Kim, N.-S. Kim, H. Lee, “Deflection
estimation of a full scale prestressed concrete girderusing long-gauge fiber optic sensors.” Construction and
Building Materials, 22 (2006) 394-401.
[2] S. Vurpillot, G. Kreuger, D. Benouaich, D. Clement,
D. Inaudi, “Vertical deflection of a pre-stressed concrete
bridges obtained using deformation sensors andinclinometer measurements.” ACI Structural Journal,(1998) 518–526.
[3] J.R. Casas, P.J.S. Cruz, “Fiber optic sensors forbridge monitoring.” Journal of Bridge Engineering, 8:6(2003) 362–373.

― 85 ―

3D STRUCTURAL ANALYSIS MODEL AND EARTHQUAKE RESPSONSE ANALYSIS OFWOODEN FIVE-STORIED PAGODA
I.TSUWA1, T.NAKAGAWA2, and M.KOSHIHARA3
1

Structural Engineer, The Japanese Association for Conservation of Architectural Monuments, Tokyo, Japan,
2
Senior Researcher, National Institute for Land and Infrastructure Management, Ibaraki, Japan,
3
Professor, Institute of Industrial Science, the University of Tokyo, Tokyo, Japan

Keywords: traditional timber structure, five-storied pagoda, non-linear dynamic behavior, Kumimono, seismic
performance, rocking motion,

2. ANALYSIS METHO
We developed a three dimensional analysis model like
Figure 1 and analyzed with the analysis soft, wallstat,developed by BRI in Japan. An analysis model was 1/5
scale model of wooden five-storied pagoda used forshaking table tests in an earlier research [2]. The
evaluation method of vertical stiffness of Kumimono wasdeveloped based on the compression tests. We applied themethod to the evaluation of roof truss. Analysis results
were compared with shaking table tests ones.
3. RESULTS
The analysis results almost agreed with experimental oneson the time history of horizontal displacement in eachfloor shown in Figure 2. However there are somedifferences in the largest deformation upper the 4th floor
between the experiment and analysis. It was seen that the
uplift of a plate under columns upper the 4th floor was
small in the analysis result compared to the experimental
one.

1: Rotational stiffness of restoring force of columns
2: Rotational stiffness of column-Nuki joints
3: Shear stiffness of wall
4: Shear stiffness of Kumimono
5: Vertical stiffness of truss roof
6: Vertical stiffness of Kumimono

6.607m1. INTRODUCTIONS
The conservation and utilization of traditional woodenbuildings have been conducted actively in Japan. Even
existing such buildings need sufficient seismic
performances because many general people visit and use
them. Therefore the quantitative evaluation of seismic
capacity for traditional timber structures is indispensable.
However most of them are built based on the knowledge
and experience of carpenters. We need more experimental
and analytical data when we evaluate the seismic
performance of such buildings quantitatively.
This study focused on the five-storied pagoda in
traditional wooden buildings. It is said that no pagoda
have collapsed against large earthquakes and they have agood seismic performances. Studies on wooden fivestoried pagoda have been increasing. However it is notclarified how a wooden pagoda can be modeled withtheories or the formularization of restoring force
characteristics proposed previously.
An early research [1] pointed the effect of the bending
deformation of a whole building on the vibration
characteristics of a pagoda. It said that a bending
deformation was caused by the rocking behavior of eachfloor, and compressive strain inclined to the grain will
happen at the parts of Kumimono etc., when each floor
moved rocking motion. Kumimono is called bracket
complexes in English and is between columns and roof.
In this study, we focused on the vertical stiffness of
Kumimono and roof truss as the structural elements
causing rocking motion and developed the three
dimensional analysis model of a five-storied pagoda
considering rocking motion. In addition earthquake
response analysis was conducted.

6
5
3
2
1

4

1.296mFigure 1. Analysis model
Experiment

Analysis

1st floor

Top

Figure 2. The time history of horizontal displacement
4. CONCLUSIONS
The analysis model could almost estimate anexperimental result. In order to increase the accuracy ofanalysis model, we need to reconsider the setting methodof the spring under the plate furthermore and to make a lotof parametric studies.
REFERENCES[1] N. Kawai, et al., 2007. Research on Dynamic
Characteristics of Five-story Pagodas. The report ofGrants-in-Aid for Scientific Research between 2004-2006.
[2] K. Chiba, et al., 2005. Seismic performance of fivestoried pagoda based on scaled model experiment part 18., Summaries of Technical Papers of Annual Meeting
Architectural Institute of Japan, C-1, pp.481-482.

― 86 ―

CODE2015

Technical Program – Day 2

Session 10

GENERALIZED FINITE ELEMENT FORMULTATION ON FIBER BEAM ELEMENTS FOR
DISTRIBUTED PLASTICITY ANALYSIS
K. PARK1, H. T. KIM2, and D.-J. KIM3
1

Assistant Professor, Department of Civil and Environmental Engineering, Yonsei University, Seoul, Republic of Korea,
Graduate Student, Department of Civil and Environmental Engineering, Yonsei University, Seoul, Republic of Korea,
3
Assistant Professor, Department of Architectural Engineering, Kyung Hee University, Yongin, Republic of Korea,
Correspond to Prof. D.-J. KIM (djkim@khu.ac.kr)

2

Keywords: Fiber beam element, Plastic hinge analysis, Generalized finite element method, Enrichment function.
1. INTRODUCTION
In the nonlinear finite element analysis of fiber beam
members, plastic hinges start to form at the locations
where extreme fibers reach the yield strength. As loadincreases, the number of fibers that yield also increases,
thus the region of the yielded cross-sections exist over acertain length of the member. As a result, the overall
stiffness of the cross-section abruptly decreases in that
region and the deformed shape of the structure shows astrong gradient. If the standard finite element is used, thiseffect can be accurately captured only by using a finemesh, which may lead to a large computational cost.
Furthermore, the progressive collapse analysis of a frame
structure requires the modeling of plastic hinges, and, ingeneral, a large number of plastic hinges are formed
before the collapse of the structure. Consequently, thestandard finite element method is not an ideal choice toperform this kind of analysis. In this study, we propose ageneralized finite element technique that can accurately
capture the deformed shape of fiber beam members withplastic hinges. This methodology enables us to obtainaccurate solutions with a relatively small computationalcost, if compare to the standard finite element method.
2.
GENERALIZEDFINITEELEMENTFORMULTAIONFORPLASTICHINGE
MODELING
The deformed shape and its derivative of the classicalEuler-Bernoulli beam with a plastic hinge are presented
in Figure 1, as discussed in [1]. In this case, the plastic
hinge modeling is based on the distributed plasticity. The
generalized finite element method can model an internal
hinge by adding new degrees of freedom to the elementwith the plastic hinge, instead of using a more refinedmesh for the member with plastic hinges. As a result, thedisplacement of the element with a plastic hinge can be
represented by the addition of a smooth approximation(v0) and a plastic hinge approximation (vE), i.e.
v( x) v 0 ( x)  v E ( x)


2

DL

 ( x)  L ( x)u  BL ( x)u  1 i 1

i

0
i

i

E
i

(1)

where  is a partition of unity function, L i apolynomial enrichment function, and B an enrichment
function describing strong gradient in the yield region
due to the formation of the plastic hinge. Since theformulation of the Euler-Bernoulli beam requires theapproximated displacement field to retain the C1 feature,
the Hermite cubic shape functions associated withdisplacement are used as the partition of unity function in
Eq. (1). Figure 2 shows the partition of unity andenrichment functions in Eq. (1). Several numericalexamples demonstrate the effectiveness of this approach.

Figure 1. Deformed shape and its derivative of a fiber
beam member with a plastic hinge

Figure 2. Partition of unity and enrichment functions forplastic hinge modeling
3. CONCLUSION
In this study, we briefly present a generalized finite
element formulation for the modeling of fiber beam
elements with plastic hinges. Proper enrichment
functions are also proposed to accurately capture thedeformed shape of the beam with plastic hinges.
REFERENCES[1] D.-J. Kim, J. Park, C.A. Duarte, S.-G. Hong, EfficientFirst-order Plastic Hinge Analysis Based on the
Generalized Finite Element Method, 11th WCCM, (2014),
Barcelona, Spain.

― 87 ―FINITE ELEMENT FORMULATION FOR FUNCTIONALLY GRADED BEAM STRUCTURES
K. YOON1, P.S. LEE2, and D.N. KIM1,3
1

Department of Mechanical and Aerospace Engineering, Seoul National University, Gwanak-ro 1, Gwanak-gu, Seoul
08826, Republic of Korea,
2
Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology, 291 Daehak-ro,
Yuseong-gu, Daejeon 34141, Republic of Korea,
3
Institute of Advanced Machines and Design, Seoul National University, Gwanak-ro 1, Gwanak-gu, Seoul 08826,
Republic of Korea,
Correspond to Prof. D.N. KIM (dnkim@snu.ac.kr)

Keywords: Finite element method; Beam; Functionally graded material; Torsion; Warping
1. INTRODUCTION
Functionally graded materials (FGMs) have received
great interest from engineers and researchers as apotential structural material for future engineering
applications. So far, a large number of investigations
concerned with the analysis of FG beams has beenreported. However, most research work has been limited
to analysis on 2D plane.
In this study, we present a formulation of three
dimensional FG beam finite element for predicting
structural responses. Here we briefly introduce theformulation of beam element and warping effects forfunctionally graded material.
2. FORMULATIONS
The geometric interpolation of the continuum
mechanics based beam finite element including warping
displacements are written as,
t

x(m)

q

q

k 1

k 1

The governing equation of the warping function for
FG beams are given by§ w 2 f k( m ) w 2 f k( m ) · wG ( m ) § wf k( m )
(m) ·
¸
¨
¸

G ( m ) ¨¨22¸¨ wyˆ  zˆ ¸
ˆ
ˆ
ˆ
w
y
w
w
y
z
©
¹
©
¹


q

t

k 1

0

,
(2)
where G (m ) is a functionally graded shear modulus.
Applying some algebraic manipulation to Eq. (2), thefinite element formulation of the warping function for
FG beams are derived asn ªº§ wGf ( m ) wf ( m ) wGf ( m ) wf ( m ) ·¦ «³: G ( m) ¨¨ wyk wky  wzk wkz ¸¸d: ( m) »
m 1«
»¼
©
¹
¬
(m)

¦ >³* Gf k( m)G ( m) n (ym) z ( m)  nz( m) y ( m) d* ( m) @
n

m 1

(m)

n ª
wG ( m ) ( m ) ( m ) º ¦ « ³ Gf k( m )z d: »:wym 1¬¼¦ hk (r )t x k  ¦ hk (r ) yk( m) t Vyk

(m)

n ª
wG ( m ) ( m ) ( m ) º ¦ « ³ Gf k( m )y d: »
:
wz
m 1¬
¼.

q

 ¦ hk (r ) zk( m ) Vzk  ¦ hk (r ) f k( m ) tD k t Vxk ,·wG ( m ) § wf k( m )¨ yˆ ( m ) ¸¸wzˆ ¨© wzˆ
¹

(1)

k 1

in which hk (r ) is a longitudinal directional
interpolation function, t x k is a position vector for beam
node k , t Vxk , t Vyk and t Vzk are directors, yk(m ) andzk(m ) are discretized cross-sectional position function,and f k(m ) and tD k are warping function and warping
degree of freedom, respectively.

Figure 1. Concept of the continuum mechanics based
beam element.

(m)

(3)
By using the numerical warping function obtained
from Eq. (3), the accurate and reliable twisting response
can be predicted. The detailed derivation is welldescribed in Ref. [1, 2, 3, 4].
3. CONCLUSIONSIn this paper, we presented a finite element
formulation of FG beams. Since the warping function
appropriately calculated by using newly derived
equations, the element can accurately predict twisting
actions.
REFERENCESG[1] K. Yoon, P.S. Lee, D.N. Kim, Compos. Struct. 132(2015) 1231–1247.
[2] K. Yoon, Y.G. Lee, P.S. Lee, Struct. Engrg. Mech. 43(2012) 411-437.
[3] K. Yoon, P.S. Lee, Comput. Methods Appl. Mech.
Engrg. 281 (2014) 106-130.
[4] K.J. Bathe, Finite element Procedures, Prentice Hall,1996.

― 88 ―ACOUSTIC WAVE SIMULATION BASED ONGPSEUDO IMPULSE RESPONSE
T. YOSHIMACHI1, M. TANIGAWA2, and K. KASHIYAMA3
1

Graduate student, Department of Civil and Environmental Engineering, Chuo University, Tokyo, Japan
2
Senior research engineer, Center for Social System Engineering, Shimizu Corporation Tokyo, Japan
3
Professor, Department of Civil and Environmental Engineering, Chuo University, Tokyo, Japan
Correspond to Prof. K. KASHIYAMA (kaz@civil.chuo-u.ac.jp)

Keywords: Noise,Impulse response, CIP Method, AMR, Convolution
1. INTRODUCTION
The evaluation of noise is important for the planningand designing of traffic facility, construction work and soon. The acoustic wave theory is a suitable method to
obtain the accurate numerical solutions because it is easy
to treat the complicated geometry and boundary
conditions. However, the method based on the acoustic
wave theory require a lot of CPU time. The convolutionof impulse response and input sound is a suitable method
to obtain the output for various type of sound source.
This paper presents a method for acoustic wave
simulation. The CIP method using adaptive meshrefinement [1] is employed for solving the impulse
response. The present method is applied to a threedimensional wave propagation problem in order toinvestigate the validity of the method.
2. GOVERNING EQUATIONS AND NUMERICALMETHODThe equation of motion and the equation of continuity
in an acoustic medium are employed for the governing
equations. The CIP method, which is a highly accurate
method for solving advection equations, is employed forthe discretization.
The adaptive mesh refinement method is employed torealize the efficient simulation. The convolution
quadrature developed by Lubich [2] is employed for thepseudo-impulse response. The pseudo-impulse generated
by this method has smooth waveforms and almost flat
frequency characteristics.
3. NUMERICAL EXAMPLE
The present method is applied to the three-dimensional
wave propagation problem to investigate the validity andefficiency. The computational model is shown in Figure1(left). The virtual sound source is arranged at the point
which is 4 m away from the flat boundary. The receiving
points are arranged at the points which are 8 m and 16m
away from the virtual sound source.
Figure 1(right). shows the sound propagation at
500step. From this figure, it can be seen that the pseudo
impulse from the flat boundary propagate in the threedimensional space. Figure 2. shows the sound pressure at
the receiving points 1 and 2. The maximum soundpressure at receiving point 1 is 0.3514Pa, and at point 2
is 0.1757Pa. The distance attenuation of sound pressure
between the receiving points 1 and points 2 is to be
0.17530.3514=0.498, and which is good agreement with

Figure 1. Sound wave propagation

Figure 2. The sound pressure at the receiving point
the theoretical values of
propagation in free space.

spherical

sound

wave

4. CONCLUSIONS
A acoustic wave simulation method based on pseudo
impulse analysis has been presented in this paper. The
convolution quadrature method developed by Lubich has
been employed for the pseudo-impulse response. The
CIP method using adaptive mesh refinement has been
usefully employed for solving the impulse response. The
solution of present method was good agreement withtheoretical values. The verification of the present method
to the actual noise problem is left for the future work.
REFERENCES[1] Y. Moriya, M. Tanigawa and K. Kashiyama, A study
on sound wave analysis based on CIP method using
adaptive mesh refinement㸪Journal of Japan Society of
Civil Engineers, Ser. A2 (Applied Mechanics (AM))
69(2),I_195-I_ 202, 2013㸬
[2] C.G Lubich, Convolution quadrature and discretized
operational calculus. I. Numerische Mathematik, Vol. 52,
pp.129-145, 1988,― 89 ―EVALUATION OF CONTACT AREA ON KNEE JOINT CARTILAGE USING IN-VIVO SUBJECTSPECIFIC FINITE ELEMENT MODEL
K.T. KANG1, S.H. LIM1 and H.J. CHUN1
1Intelligent Structure and Integrated Design Laboratory, Department of Mechanical Engineering, Yonsei University,
Seoul, Republic of Korea,
Correspond to Prof. H.J. CHUN (hjcun@yonsei.ac.kr)

Keywords: Finite element analysis, in-vivo validation, Subject-specific knee model, Weight-bearing MRI
1. INTRODUCTION
Contact pressure increases or alterations in knee joint
structure and contact area are critical causes of knee pain
and osteoarthritis [1]. This is an active area of research,
the exact mechanical behavior of the knee joint and thecauses of knee joint injury have not been completely
elucidated. Carey et al. have developed a technique tovalidate the Finite Element(FE) model in-vivo using
dynamic stereo-radiography data [2]. However, it
requires a dual-fluoroscopy apparatus or weight-bearing
MRI. The present FE study is cost effective and efficientand allows us to resolve problems as they develop, which
is not generally the case in experimental studies.
2. MATERIAL AND METHOD
The subject was positioned supine with his right knee
at full extension in an axial loading device. Images wereacquired using a 3.0 Tesla Magnetic ResonanceImaging(MRI) scanner and a custom designed knee joint
cadio coil. A 3D FE model of a healthy lower extremity
was developed from Computed Tomography(CT) images
obtained with a light-speed volume CT scanner. An
Iterative Closet Point(ICP) algorithm was used foralignment of X-ray and supine MRI(XSM).

Figure 1. Alignment of X-ray and supine MRI using an
ICP algorithm

Figure 2. Interaction volume between non-deformed
tibiofemoral cartilage: (a) X-ray and supine MRI (XSM);(b) Weight-bearing MRI(WBM)
3. RESULTThe XSM and WBM methods were employed forverification of the subject-specific FE model. Table 1
shows the values of contact area and contact deformation
for subject specific XSM method, WBM method and FE
model.

Figure 3. Results of cartilage contact area and contact
deformation distribution: (a) X-ray and supine MRI(XSM); (b) Weight-bearing MRI(WBM); (c) FE model
Table 1. Result of cartilage contact area and contact
deformation
Contact
Contact area
deformation
(mm2)
(mm)
Medial Lateral Medial Lateral
X-ray and
224.1
198.9
0.51
0.43
supine MRI
Weight228.8
193.4
0.49
0.41
bearing MRI
FE model
246.8
188.5
0.53
0.39
4. CONCLUSIONS
Verification was conducted by comparing the resultsfrom the model with either literature data or in vitro
cadaveric data. Our study reduced the error in six
degrees of freedom in alignments of 2D and 3D images.
The errors between the two images were 0.043 ± 0.011
and 0.021 ± 0.005 in coronal and sagittal views,
respectively. The XSM method used in this study is anacceptable validation method in general situations andshowed very similar results to the WBM validation
method for the in vivo model. In conclusion, the XSM
method proposed in this study represents a step forward
in subject-specific in vivo validation of original subjects
in biomechanics.
REFERENCES[1] J.J. Elias, D.R. Wilson, R. Adamson, A.J. Cosgarea, J.
Biomech. 37 (2004) 295-302.
[2] R.E. Carey, L. Zheng, A.K. Aiyangar, C.D. Harner, X.
Zhang, J. Biomech. 136 (2014) Paper No.041004.

― 90 ―NONLINEAR FLEXIBLE MULTIBODY APPROACHES USING CO-ROTATIONAL FINITEELEMENTS
H. CHO1, H.S. JOO1, C.W. PARK1, and S.J. SHIN2
1

Graduate student, Department of Mechanical and Aerospace Engineering, Seoul National University, Seoul, Korea,
2
Professor, Department of Mechanical and Aerospace Engineering, Seoul National University, Seoul, Korea,
Correspond to H. CHO (nicejjo@snu.ac.kr)

Keywords: Flexible multibody analysis, Co-rotational finite elements, Lagrange multiplier
1. INTRODUCTION
Flexible multibody analysis is generally required withgrowing demands in the mechanical engineering field. It
combines concepts and solution techniques from both
flexible structural mechanics and the rigid bodymechanics. During the last few decades, the researchers
established existence and uniqueness of solutions andbenign numerical methods related with the multibody
analyses. The Lagrange multiplier technique plays a majorrole in modeling the multibody analysis with constraints[1]. In this paper, an improved and more versatile
structural analysis applicable for nonlinear flexible
multibody system is developed. In more detail, nonlinear
finite element based on co-rotational (CR) framework isdeveloped. Then, interconnecting methods using the
Lagrange multipliers are implemented to constrain the
kinematic relationship among the flexible components.
Moreover, parallelized algorithm of the analysis will besuggested.
2. CO-ROTATIONAL FINITE ELEMENTS
The CR framework is independent of the specific local
element formulation. Thus, a number of existing robust
and accurate geometrically linear plane elements can bereused and extended to geometrically non-linear analyses.
A basic idea of this approach is to extract pure elastic
deformation from the total motion using the CR frame
which translates and rotates along the element.
Accordingly, three kinds of coordinates are to be defined
to deal with such decomposition of the behavior. By
following the transformation between the coordinates, theglobal stiffness matrix and internal force vector of theelement can be defined. Details of its derivations arepresented in Ref. [2].
3. MULTIBODY APPROCHESIn this section, both serial and parallel versions ofmultibody analysis are presented. Classical Lagrange
multiplier method is employed to enforce the kinematicconstraint among the components. This method yields astraightforward formulation, in which addition of thevirtual work of individual bodies is performed to obtainthe total virtual work of entire bodies, considering each of
the kinematic constraints between any two bodies or two
adjacent nodes by the Lagrange multiplier vector. To
ensure computational robustness and efficiency, theparallelized algorithm based on the improved finite
element domain decomposition method using local
Lagrange multipliers is also developed [3]. And it isapplied for the present analysis. Currently, static algorithm
regarding nonlinear analysis is developed. Thus, it will be
extended to nonlinear transient analysis.
4. NUMERICAL RESULTS
To examine the present analysis, interconnected threebeam structure under a tip load is considered (Fig. 1(a)).
The CR planar element is employed and both serial andparallel analyses are conducted. The deflectedconfiguration is illustrated in Fig. 1(b). To verify anefficiency of the present parallel analysis, scalability test
is conducted (Fig. 2).

(a) Configuration of multibody(b) Deformed
analysisconfiguration
Fig 1. Interconnected beam structures under tip load

Fig 2. Scalability test of the present parallel analysis
5. CONCLUSIONSIn this paper, flexible multibody analysis using the CR
elements is developed. Extension of the present analysis
in order for transient analysis will be conducted.
REFERENCES[1] B. Simeon, On Lagrange Multipliers In Flexible Multibody Dynamics, Comput. Methods Appl. Mech. Eng. 195(2006).
[2] C. Felippa and B. Haugen, A Unified Formulation of
Small-strain Corotational Finite Elements: I. Theory,
Comput. Methods Appl. Mech. Eng. 194 (2005) 22852335.
[3] J.Y. Kwak, T.Y. Chun, S.J. Shin, O.A. Bauchau,
Domain Decomposition Approach to Flexible Multibody
Dynamics Simulation, Comput. Mech. 53 (2014) 147-158.

― 91 ―

APPLICATION OF TRANSFORMED AUTOCORRELATION FOR STRUCTURAL SYSTEMSWITH LOGNORMALLY DISTRIBUTED RANDOM VARIABLES
H.C. NOH1
1

Associate Professor, Dept. of Civil and Environmental Engineering, Sejong University, Seoul, S. Korea,
Correspond to Prof. H.C. Noh (cpebach@sejong.ac.kr)

Keywords: Response variability, Lognormal random variable, Weighted integral method, Uncertainty, Correlation
1. STRUCTURAL UNCERTAINTIESThe uncertainties or randomness in structural systems
are caused by a variety of reasons. Probabilistic
assumption on the systems is the first source of
uncertainties from the perspective of uncertain responses.
However, this source is intrinsically unpredictable and isbound to have errors. The uncertainty which is natural,
influential and conceptually easy to implement is thematerial randomness spread out over the domain of thesystem. The other one is categorized as geometricalparameter which characterize the shape of the system.
In order to deal with the randomness in the structural
systems, we have two types of methods. One is the
MCS(Monte Carlo Simulation) which can be called astatistical method and the other is indirect method whichuses indirect function of autocorrelation which implies thestatistical characteristics of the given random fields. The
MCS is easy to implement if we see this method from theviewpoint of analysis itself since the essence of thismethod is repetitive application of finite element analysis.
MCS, however, is based on the actual (theoretically)
random fields which are numerically generated relying on
the various random field generation schemes. We can use,for example, spectral representation and/or statisticalpreconditioning etc. In the case of indirect method, it hasa relative advantage to the MCS in the computation timesince it deals only with the autocorrelation (or spectral
density) functions which represent the random field in an
indirect way. The further advantage of the indirect method
to the MCS is that we can obtain a regular results which is
not random. With this, we mean the results of MCS isirregular depending on the condition of analysis.
Whatever the scheme of analyses, almost all the randomparameters have been assumed to be distributed in a
Gaussian way [1]. (Not all, but almost.) It is natural tothink that a random parameter is to have probability
density of Gaussian, and the central limit theorem isgiving us legitimacy for this assumption. The several
practical random parameters, however, show probability
density other than Gaussian, e.g. lognormal distribution,
beta distribution etc. To assess the uncertain response due
to random parameters other than Gaussian, the random
fields need to be numerically generated. For this purpose,the so-called translation theory can be employed, withwhich the Gaussian random fields (and/or processes) canbe transformed into the fields distributed in the other way.
Accordingly, we need to follow the following procedure:
the generation of Gaussian fields, and then transform theminto the target fields distributed in a nonGaussian way, andthen perform N repetitive analyses.

In the case of indirect method, however, direct
application of nonGaussian random parameter is noteasily achieved since the autocorrelation (or spectral
density) function is not easy to get. Fortunately, we canderive the corresponding autocorrelation function for the
lognormally distributed random parameters.
In this study, by using the autocorrelation function
transformed from that for Gaussian random variable isadopted to apply to the stochastic weighted integral
formulation.
2. FORMULATION
The form of autocorrelation function for lognormal
variable is as follows [2]:
CLN   LN exp CN   1(1)
Where, three parameters of  ,  ,  are determined



 

based on the statistics of lognormal variable, andCN denotes the correlation structure of Gaussianvariable.
If we employ Eq. (1) to be the correlation structure
for the lognormal variable, we can get the weighted
integral formulation as follows:
Ne

Ne

e1

e2

Cov  U   K LN1  E Ce1e 2 K LNT

(2)

2
E Ce1e 2      LN
CLN    k   UUT k   dV2 dV1 (3)
e1

e2 T

V V
Finally, with
Eqs. (1-3) we can get the variabilityof the response for the given system [1].
1

2

3. CONCLUSIONSHow to evaluate the response variability of structural
systems with lognormal random parameters is addressed.
By comparison with the results obtained for the Gaussian
random variable, we can get the characteristics of practicalvariability caused by the lognormal random variable.
Furthermore, by employing the translation theory, thismethod is expected to have more broad application to thepractical structures.
REFERENCES[1] Hyuk-Chun Noh and Taehyo Park, Computer Methods
in Applied Mechanics and Engineering, Vol. 200, Issues
29-32 (2011), 2397-2406.
[2] Gasper Z., Andrej T., Donald L. S. and Roberto C.,
Nuclear Instruments and Methods in Physics Research A,
727 (2013), 33∼39.

― 92 ―