Schedule for Presentation Meeting of Master-Degree Thesis in 2020
International Course of Maritime & Urban Engineering
Date :
Venue :
No.

Time

Name

1

09:30 – 09:55

AHMED Murshalin

Sanada2

09:55 – 10:20

NER Nikko Torres

Kita

3

10:20 – 10:45

RAHMAN Tahmina

Kita

4

10:50 – 11:15

5

11:15 – 11:40

NALLARAJAH
Renee Josiah
BUERLIHAN
Ayang

Supervisor

Title of thesis
FE Modelling of Reinforced Concrete Exterior Beam-Column Joint with Deficient BeamLongitudinal Reinforcement Anchorage
Flood Coping Mechanisms in In-city Resettlement Areas:
A Case Study of Barangay Malanday in Marikina City, Philippines
A Study on the Potentials of Green Blue Infrastructure (GBI) Network for Dhaka City in
Bangladesh

ChairmanMatsubara
Matsubara
Matsubara

Araki

Impacts of Tsunami Driven Debris on Onshore Coastal Structures

Araki

Tsutsumi

Study on the effect of weld bead shape on fatigue strength of Out-of-Plane Gusset

Araki

6

13:00 – 13:25

SUBRAMANIAM
Sreenath7

13:25 – 13:50

ISNAINI Rodhiatul

Iijima8

13:50 – 14:15

Thein Lin AungMa

9

14:15 – 14:40

RACHMAN Dimas
MaulanaMaki

14:55 –

at S1-312

Meeting for Evaluation:

July 29th (Wednesday) starting from 09:30
S1-312 Lecture Room

Umeda

Lunch Break
Systematic Study of Water-on-Deck Effect of Offshore Supply Vessel Running in Stern
Quartering Waves
Still Water Bending Moment Evaluation on Beam through Cargo Loading Prediction using
Kalman Filter
Development of Iso-Geometric Analysis Software JWRIAN-IGA for Weld Bead Geometry
Design Using Bayesian Optimization
Study of Automatic Berthing by Evolution Strategy with Covariance Matrix Adaptation(CMA-ES) in Real Time

Meeting for Evaluating Master-thesis Presentations, and Steering Committee Meeting

From 14:55 at S1-312 Lecture Room
Members are Steering committee members and the supervisors

Maki
Maki
Maki
Maki

International Program of Maritime & Urban Engineering / Abstract of Master Thesis

July, 2020

FE Modelling of Reinforced Concrete Exterior Beam-Column Joint
with Deficient Beam Longitudinal Reinforcement Anchorage
Murshalin Ahmed
Concrete Structure Laboratory, Department of Architectural Engineering
Key Words: Bond-slip, Deficient anchorage, Exterior joint, FEM analysis, Low strength concrete

Concretecompressive
strength/Steelyield stress
Young’s
Modulus
Axial force

where, 𝜏𝑢 is bond strength; 𝜎0 is axial stress of column; 𝜎𝐵 isthe compressive strength of concrete.
Bond slip corresponding to the bond strength was assumed as 1.0mm, referring to the study by Elmorsi et al. 5).
5-D13
10-D13Concrete11.0 N/mm2

D13342 N/mm2

D6408 N/mm2

Concrete
D13D6

8.7 kN/mm2
174 kN/mm2
171 kN/mm2
72 kN

2.2 Bond-Slip Characteristics
Naganumas' model2) was used for defining the bond stress-slip
characteristics between rebar and concrete. Two models wereconsidered for peak bond strength. In model 1 peak bond strength
in the beam and the column sections were calculated according to
the Japanese Design guidelines3) (shown in Eq. 1). For model 2,

D6 @ 60

D6 @ 60

3-D13
210
Beam section

270
Column section

2.1 FE Modelling
The specimen analyzed was the specimen J1 from the previous
study1), as shown in Fig. 1. Table 1 shows the specification of
specimen J1. An image of analytical modeling is shown in Fig. 2.
Concrete and pin supports were modeled by four-node isoparametric quadrilateral elements in 2-D discretization.
Longitudinal bars in the beam and column were modeled bydiscrete truss elements to consider bond behavior between
reinforcing bar and concrete. The bond behavior was modeled byintroducing a linkage element on the interface between the truss
element and the concrete element. All other reinforcing bars weremodeled by embedded elements, assuming perfect bonding. The
mechanical properties of concrete and reinforcements werereferred to the experimental results, as shown in Table 1.
Table 1 Specification of the specimen
bxD
270 mm x 270 mm
ColumnMain rebar
10-D13
Shear reinf.
D6@60
bxD
210 mm x 370 mm
Beam
Main rebar
8-D13
Shear reinf.
D6@60

(1)

𝜎𝐵

Straight Anchor

235

2100

2. Finite Element Analysis𝜎𝜏𝑢 = 0.7 (1 + 0 ) 𝜎𝐵 2/3

370

The application of straight anchorage at the exterior beamcolumn joint is often found in developing countries because of thesimplicity of its installation. In this study, nonlinear 2-D finite
element analysis was conducted on reinforced concrete (RC)exterior beam-column joint specimen with deficient beam rebar
anchorage and the seismic capacity was predicted. The results by
the finite element analysis were also compared with theexperimental results. The research describes the effectiveness ofmodeling of the straight anchorage considering the bond-slip
relationship between the rebar and concrete and its limitations.

simple pullout test was performed on pre-installed rebar in concrete
blocks with variable embedment length4). The peak bond strength
value found from that test was considered for the model.

270

1. Introduction

Strain gauge
unit: mm
1125

Figure 1 specimen J1 details
Axial Force
Cyclic force
-

+

5@52

(Rigid)

5@41

10@40

(Rigid)
41
35
41
4@34
41
35
41

Roller

2@46
6@55

6@50

2@48

2@42

2@44

Reinforcement(as truss element)(Rigid)Pin
41

4@47 41

Figure 2 FEM element division and boundary conditions

3. Results and Discussion
Experimental results1) in Fig. 3 showed that specimen J1 suffered
pullout at -2% drift ratio in negative loading direction. However,on the positive loading direction joint shear failure was observed.

direction, which was not observed during the experiment1). The
bond strength applied to the analysis for model 1 followed theJapanese Design guidelines3), which might be conservative
compared to that in the experimental specimen. The bond strength
applied to the analysis for model 2 is from pure pullout test, whichmight be radical compared to that in the experimental specimen.

Story Shear Force, Q (kN)

30
Experiment

20
10
0

3000

-10
-20
-3

-2

-1
0
1
2
Story Drift Ratio, R (%)

3

4

Figure 3 Experimental results1)
The story shear force vs. story drift angle relationships from the
FE analysis and experimental study are compared in Fig. 4. Model1 reached a maximum strength of +25.5 kN & -19.0 kN in thepositive and negative directions, respectively. Whereas for the caseof model 2 it was +27.6 kN and -21.5 kN, respectively. A large
slip/pullout was observed in both positive and negative loading
directions for model 1. In the positive loading direction beamrebars slipped at R = +2.0% drift cycle. While, in the negative
loading direction beam rebars slipped at R = -1.5%, as shown in
Fig. 5. No rebars yielded in the FEM analysis and strains observedin the corresponding cycles were lower than those of theexperiment. On the other hand, pullout was not observed in either
loading direction for model 2. Rebars yielded both in positive andnegative direction at 1.5% drift ratio, as shown in Fig. 5.

Strain (µ)

2000

-4

Rebar Yielding

1000

0

-1000

Experiment
FEM model 1
FEM model 2

-4

-3

-2

-1
0
1
2
Story Drift Ratio, R (%)

3

4

Figure 5 Beam-column joint end rebar strains
Figure 6 shows the damage to the specimen of model 1 at R= 1.5% drift ratio when pullout was observed along with the damageof the experimental specimen J1.

40

Story Shear Force, Q (kN)

30

Experiment
FEM model 1

Beam Yielding(a)(b)
Figure 6 Damage at R= -1.5%; (a) experiment, (b) FE analysis

20
10

4. Conclusions

0
-10
-20
Beam Yielding

-30

-4

-3

-2

-1
0
1
2
Story Drift Ratio, R (%)

3

4

40

Story Shear Force, Q (kN)

30

Experiment
FEM model 2

Beam Yielding

20
10
0

References-10-20
Beam Yielding-30

The following findings were obtained from this study:
1. Finite Element analysis showed pullout failure of beam
longitudinal rebar with straight anchorage can be simulated
by considering bond behavior between reinforcing bar andconcrete.
2. Model 1 underestimated the joint capacity when compared tothe experiment because the peak bond strength used wasfrom design guideline, which provides conservative
estimations.
3. Model 2 overestimated the joint capacity when compared tothe experiment because the peak bond strength used wasfrom element test. In real condition, the joint is damaged andthe acting bond strength is lower than the simple element test.
4. To get a better estimation of the joint characteristics, bondstrength needs to be determined form the sub-assembly test.

-4

-3

-2

-1
0
1
2
Story Drift Ratio, R (%)

3

4

Figure 4 Story shear vs. story drift ratio relationships
For model 1 the FE analysis was successful in predicting thepullout failure mode in the negative loading direction. Furthermore,
the FE analysis also showed pullout in the positive loading

1) Wardi, S., Saha, N., Sanada, Y., and Takahashi, S. (2018),
Proceedings of AIJ Annual Convention 2018, Paper No. 23334.
2) Naganuma, K., Yonezawa, K., Kurimoto, O. and Eto, H.:, 13th
WCEE, Paper No.586, August., 2004.
3) Architectural Institute of Japan: Design Guidelines for
Earthquake Resistant Reinforced Concrete Buildings Based on
Inelastic Displacement Concept (Draft) 1997.
4) Saha, N., Master Thesis, Osaka University, 2019.
5) Elmorsi M, Kianoush MR, and Tso WK., Can. J. of Civil Eng.
2000; 27: 490–50

International Program of Maritime & Urban Engineering / Abstract of Master Thesis

July 22, 2020

Flood Coping Mechanisms in In-city Resettlement Areas:
A Case Study of Barangay Malanday in Marikina City, Philippines
Nikko Torres Ner
Architectural and Urban Planning Lab, Department of Architectural Engineering
Key Words: in-city resettlement, coping mechanisms, urban flooding, vulnerability, disaster risk management

1. Introduction
A World Wide Fund for Nature (WWF) study in 2009 rated
Metro Manila as the most at risk among 11 mega-cities due to its
exposure to typhoons and flooding (Shaw 2010). This problem iscompounded by rapid urbanization and decades of unplanneddevelopment which had pushed the marginalized sector to live inblighted and high-risk areas in Metro Manila.
A 2017 study found that in-city resettlement was more efficientin the use of government resources in addressing the housing
challenges in the long term compared to off-city resettlement(Ballesteros 2017). Marikina City has been implementing an
in-city resettlement program since 1992. The study area,
Barangay Malanday, is one of the six relocation areas for
Marikina City’s in-city resettlement program (See Fig. 1). The
program was initiated with the motivation of revitalizing the
Marikina River and relocating informal settlers away fromidentified danger zones at the time (Iglesias 2008).
This study intends to examine Barangay Malanday, an in-cityresettlement area at risk from the effects of urban flooding, andinvestigate its flood vulnerability and coping mechanisms to build
on existing qualitative studies that assess local knowledge andbest practices towards safer cities.

2.2 Methodology
A survey which covered 50 residents was conducted in March
2020 to gather data about respondents’ area’s socioeconomic
profile and investigate the various flood coping mechanisms
practiced by the residents. Since the area experiences flooding
frequently, questions were limited to experiences pertaining to the2009 and 2018 floods. Additional key informant interviews wereconducted with homeowners’ association leaders about theircommunity’s experiences in the 2009 and 2018 floods. Throughsemi-structured interviews, insights into their coping strategies as
a community were recorded and transcribed.
Map analysis was used to gain a better understanding of thestudy area’s level of exposure to flooding. Mapping activities
were also performed to observe the effects of the respondents’
coping strategies on the damage to their real and personal
properties.

3. Results and Discussion
3.1 Exposure
3.1.1 Flood Hazard
The study area is situated in low-lying elevation (8-12 masl)
and is in close proximity to the Marikina river. As a result, it ishighly susceptible to riverine flooding caused by heavy rainfall
especially in the monsoon season (June to August). The entiretyof the study area falls under the very high flood susceptibility
category, or areas likely to experience flood heights of greaterthan two meters as seen in Figure 2.

Fig. 1 In-city resettlement area locations (Iglesias 2008)

2. Research Design
2.1 Framework of Analysis
The investigation into the study area’s vulnerability and coping
mechanisms was guided by the Intergovernmental Panel on
Climate Change’s (IPCC) framework of vulnerability:
Vulnerability= (Exposure + Sensitivity)-Adaptive Capacity• Exposure is the nature and degree to which a system
is exposed to significant climate variations.
• Sensitivity is the degree to which a system is affectedby climate variability.
• Adaptive Capacity is the ability to adjust to climatechange to moderate potential damages

Fig. 2 Flood hazard map of the study area
3.1.2 Climate Change Projections
Models project an average temperature increase of 1.8 to 2.1
Degrees Celsius in 2050 (PAGASA 2011). With this projected
temperature increase, it is expected that storms will also increase
in intensity and cause more deleterious effects on the mostvulnerable populations. The amount of rainfall in Marikina City
during the monsoon season (June to August) is also expected toincrease by 21.3% (PAGASA 2011) from 1971-2000 levels andwill increase the probability of intense flooding in the study areain the future.
3.2 Sensitivity
Three economic classes in the study area were identified fromthe survey, namely: poor, low-income, and lower middle class.
Majority of the households fall under the poor category with 23
respondents saying they have a monthly household income of lessthan Php10,000 followed by low-income (17) and lower middle
class (10) (See Table 1). The classification is based on criteria set
by the Philippine Institute for Development Studies (PIDS) whichlooks at per capita income to assess standard of living.
Majority (46%) of the respondents said their household size
was composed of 4-6 persons. This is consistent with the average
household size in Marikina City which is 4.59 in 2019.
Table 1 Socio-economic profile (N=50)
Household Size
Count
1-3
9
4-6
23
7-9
14
<10
4
Monthly Household Income (PhP)Count
Poor (>10,000)
23
Low-income (10,001-20,000)
17
Lower Middle Class (20,001-40,000)
10

%
18
46
28
8
%
46
34
20

3.3 Adaptive Capacity
3.3.1 Housing Modifications
Respondents were asked what modifications, if any, had they
done to their houses to protect against flooding following the2009 flood. They were then asked to rate and describe the degreeof damage to their houses during the 2009 flood and then the2018 flood. Three types of housing modification were observed,
namely: building upper floors, rebuilding/reinforcing the housewith more durable materials, raising/strengthening of groundfloor. Figures 3 and 4 show the degree of damage to real property
to the flood in 2009 compared to the flood in 2018.

Fig. 4 Degree of damage to real property in 2018 flood
3.3.2 Behavioral Changes and Preparedness
All respondents mentioned proactively heeding early warning
systems and evacuation measures since 2009. In the 2009 flood,48 (96%) of the respondents reported suffering high damage totheir personal property because of complacency and ignoring theevacuation warnings. This heightened alertness from the residents
during the 2018 floods translated to only 4 of the 50 respondents
suffering high damage to their personal properties (See Table 2).
Homeowner’s association leaders also reported that aside fromtheir original function as facilitators of the mortgage program,they have taken on new roles such as organizing disaster training
workshops and assigning specific responsibilities to their officers
in times of disasters.
Table 2 Damage to personal property in 2009 vs 2018 (N=50)
Degree of Damage
2009 Flood2018 Flood
Low
031
Medium
2
15
High
48
4

4.

Conclusion

Significant reduction on the impact of floods on the study area
was observed in the comparison of the aftermaths of the 2009 andthe 2018 floods. Reduction in damage to personal property saw
the most positive improvement. This can be attributed to thedisaster preparedness of individual households and theircommunity’s disaster risk management initiatives. Reduction indamage to real property, however, could still see some areas for
improvement. Due to their economic situation, majority of thehouseholds can only afford incremental and low-quality
improvements to their houses. A multisectoral approach which
considers the residents’ realities and experiences could prove tobe beneficial in uplifting their current situation.

References

Fig. 3 Degree of damage to real property in 2009 flood1) Shaw, R.: CDRI (Kyoto, Japan), pp. 3-4, 2010.
2) Ballesteros, M.: Evaluating the Benefits and Costs ofResettlement Projects-A Case Study in the Philippines,
IDEAS (Exeter, United Kingdom), pp.257-271, 2017.
3) Iglesias, G.R: ADPC (Bangkok, Thailand), pp.5-7, 2008.
4) McCarthy, J. et al.: IPCC (Cambridge, UK), pp. 21-22, 2001.
5) PAGASA: Climate Change in the Philippines, PAGASA(Quezon City, Philippines), p.35, 2011.

International Program of Maritime & Urban Engineering / Abstract of Master Thesis

July, 2020

A Study on the Potentials of Green Blue Infrastructure (GBI) Network
for Dhaka City in Bangladesh
Tahmina Rahman
Laboratory of Urban planning and Design (Lab. 3),
Department of Architectural Engineering
Key Words: Green-Blue Infrastructure, ecological services, neighborhood planning, Dhaka

1. Introduction
Cities currently possess the highest proportion of the world’s
population, with an estimated 54% of the global populationt1).
With this rapid growth of the urban population, the regenerationand conservation of urban ecosystems are becoming a serious
challenge, especially in developing countries. In view of thecurrent pace of urbanization, per capita urban green spaces have
reduced tremendously in the cities, especially in the Global South.
However, Green-Blue Infrastructure (GBI), as key components of
the urban ecology, can play a vital role in providing theecosystem services. The city of Dhaka has also seen its green andblue areas gradually encroached upon and destroyed2), a situation
which has steadily deteriorated since Bangladesh gained
independence. In Dhaka, both public and private sectors arecontinuously developing natural lands within and on the outskirtsof the city, without considering the environmental impact of thisactivity by filling the water bodies or graving the greens4).
This study, therefore, sought to analyze the potential of GBI
for the urban ecological conservation and regeneration of Dhaka
city from both a macro and micro level perspective. At the macrolevel, spatiotemporal changes in green and blue infrastructures topropose a scheme for the protection and regeneration of theexisting green spaces and water bodies is analyzed. For
micro-level user’s preference to determine possible options to theadoption of new GBI in dense urban neighborhoods through adesign exercise for sustainability, connectivity, and well-being isassessed.

2. Methodology
2.1 Methods
This research adopted the case study strategy. In this study,
multiple data from different sources (Historical maps, Policy
regarding city development, Surveyed data etc.) are collected totriangulate the validity of the available information to develop an
in-depth understanding on the possibilities to develop a GreenBlue Infrastructure (GBI) Network for Dhaka City. The study was
conducted in two phases: macro-level study and micro-level study.
For the macro-level study, data was gathered through historical
map regression, policy content analysis and key informant
interviews. For the macro-level study, data was collected
through a structured questionnaire survey of 115 samples, transect
walks and field observations and photographic documentation.
The fieldwork was conducted during August and September of
2019. SPSS version 23 and Microsoft Excel were used to analyze
the questionnaire survey data.
2.2 Context
The study was undertaken in Dhaka city, the capital of

Bangladesh. Two neighborhoods were selected based on physicalmorphology (street pattern, planning organizations etc.),
availability of GBI and historical evidence of GBI integration.

Fig. 1 Dhaka Metropolitan Area and Studied Areas
Site A is a well-organized planned residential area with agrid-iron street pattern (Ward 49) with several GBI options suchas lake, lakeside park, playgrounds of different scales, roadverges, road islands and roadside green areas. Site B is anorganically developed residential area with a spontaneously
growing street pattern (Ward 40 (partial) and Ward 51) withseveral GBI options such as ponds, park, playgrounds of different
scales, road verges, road islands and roadside green areas.

3. Results
3.1 Changes in Urban Green and Blue Spaces (Macro Level)
From the historical map regression, historical dynamics ofurban waters and green spaces and traces in the connectivity ofurban waters `were analyzed.
From the analysis, a ‘large positive linear association’ wasfound between the total area of Dhaka city and total built areafrom 1924 to 2015 with correlation value r= 0.966 and R2= 0.934.
However, between the total area of Dhaka city and urban greens,the association was ‘small positive linear association’ with thecorrelation value r= 0.723 and R-squared= 0.525. and similar‘small positive linear association’ was also found between the
total area of Dhaka city and urban waters with the correlation
value r= 0.811 and R-squared= 0.658. Additionally, it was
identified that the trendline of the built area has a steep
incremental growth (Fig. 2). However, in the case of the urban
greens and urban water, a dramatic drop was observed after 1995.
Between 1960 to 1995, the proportion of the area of urban waters
was steady, but also started a sharp decline after 1995. Before
independence, the GBI area ratio with the built area was 51:49,
persisting until 1995. However, after 1995, this dramatic fallcontinued and by 2015, the loss was about 25% (from 10916 ha to8328 ha). Currently, the stock of GBI is only 27% (around 8328

ha), including restricted green and privately-owned informal
greens, mostly inaccessible to the general public.

the GBI sites was the main problem experienced by the majority(62%). Along with this, 59% of them in Site B, indicate theproblem of accessibility (either no GBI in close proximity or
nearby GBI has restricted Access). Last of all, their perception of
the benefits of having GBI in the neighborhood has been
investigated through a five-point Likert scale. A positive
perception regarding the benefits of GBI has been observed
except social benefits; 71% responded positively forenvironmental benefits, 59% responded positively for economic
benefits, 60% responded positively for wellness benefits and 48%
responded positively for social benefits.

4. Design Proposal

Fig. 2 Dhaka Metropolitan Area and Studied Areas
Moreover, building on the identified geomorphology andspatiotemporal changes from this research and comparing withthe streamlines and water flow direction of Dhaka from theprevious studies5,6), seven potential locations where the newconnections can be developed or the existing water network canbe resurfaced were identified.

In the design exercise, a projected design development has
been proposed (Fig. 4), considering Dhaka Metropolitan Building
Construction Act (2008); the rules of FAR, Mandatory OpenSpace, Maximum Ground Coverage has been applied todemonstrate how to preserve the existing GBI and introduce
interconnected green-blue spaces.

3.2 Availability of Urban and Green Spaces (Micro Level)
From the comparison of two studied areas (Fig. 3), it isevident that Site B, is denser in terms of built area (93.3%) andbelow the standard regarding the availability of GBI, (available
GBI is 6.7%, only around 10 ha). On the other hand, Site A has
27.8% of greens spaces and urban waters (around 104 ha, which
is about double of the maximum required open space of 57 ha, as
per Open Space Standards by Detailed Area Plan, Dhaka (2010)).
Therefore, Site B selected this area to demonstrate the design
proposal.
Fig. 4 Proposed Plan for Neighborhood X (NX)

5. Recommendation and Conclusion
In summary, the study showed that there has been a hugedecline in blue and green spaces in Dhaka, particularly since 1995.
However, the majority of the surveyed residents reported
immense environmental, economic and wellness benefits to theadoption of GBI in their neighborhood. Therefore, in addition to
the design proposal, the study recommends the integration ofconservation and revitalization of GBI in the urban design andplanning policy of Dhaka city.

References

Fig. 3 Surveyed Study Areas with Their GBI Options
Results from users of GBI found that the majority of the
respondents (70%) lived in apartments and 60% of the total
respondents used GBI for different activities. Among these
activities, physical exercise and social gathering were found to be
primary activities. Close proximity, multi-functionality andwell-equipped space were the main attributes of their preference.
‘General maintenance—littering, garbage and waste dumping’ in

1) Goldenberg, R., Kalantari, Z. & Destouni, G: J of LDD,
Vol.29, pp. 3607-3616, 2018.
2) Chowdhury, A. & Faruqui, S.: Physical Growth of Dhaka city.
Dhaka: Past, Present and Future, pp. 43-61, 1989.
3) Byomkesh, T., Nakagoshi, N. & Dewan, A. M. J of LEE,
Vol.8, pp. 45-58, 2012.
4) Islam, N. Dhaka Now: Contemporary Urban Development,
Dhaka, Bangladesh, Bangladesh Geographical Society, 2005.
5) Mandal, D., MT, West Virginia University, 2019.
6) Peeters, S. & Shannon, K.: Nakhara—J of EDP, Vol.7, pp.
25-46, 2011.

International Program of Maritime & Urban Engineering / Abstract of Master Thesis

July, 2020

Impacts of Tsunami Driven Debris on Onshore Coastal Structures
N. R. Josiah
Coastal Engineering Laboratory, Department of Civil Engineering
Keywords: Tsunami, wave, debris, impact force

1. Introduction
Extreme coastal disasters like tsunami always bring great
devastation to the properties and lives as witnessed in recent past
especially during 2004 Indian Ocean tsunami and 2011 Great East
Japan tsunami. Post field surveys of tsunami damages indicatedthat in addition to wave induced forces, floating objects (i.e.
debris) carried by tsunami also contributed significant amount offorces which leads to either damage or failure of the structural
facilities at inundated area (Murata et al., 2016). Debris can be ausually available objects in coastal and harbour areas can beranged from boats, shipping vessels, containers, vehicles,
damaged pieces of buildings and lumbers (Nistor et al., 2008).
Even though few design guidelines (FEMA55, 2011; FEMA P646,
2012) have discussed a need of estimation and basic approaches
to calculate debris caused impacts, a proper way to encounter thisphenomena is currently under development. This study presents
an experimental investigation conducted to examine the impact
forces exerted by different type of debris to a typical one storey
structure available in coastal area. Objective of the study is to
identify the variation of impact forces in relation to wave forces,
collision patterns and to diagnose governing parameters ofdebris-structure interaction provoked by a single debris element.

2. Methodology
2.1 Experimental facility
Wave tank facility at the Department of Civil Engineering,Osaka University was utilized to conduct the experiment. Figure
1 shows the experimental set-up prepared. Tsunami like flow was
generated as suddenly lifting the gate of the reservoir to allow thebore to propagate. Water level difference between flume andreservoir (Δh) was varied to generate different wave bores.

in Table 1. Velocity and inundation depth were measured at thelocations of V and W4 as shown in Figure 1. Velocities andinundation depth in front of the house model were in the range of
0.55 m/s to 1.02 m/s and 0.033 m to 0.074 m accordingly.
Properties of debris used in this study is presented in Table 2.
Table 1: Flow conditions
Case No.

∆h (m)

Max. velocity
(m/s)

I
II
III
IV
V
VI

0.1
0.12
0.15
0.17
0.2
0.22

0.55
0.63
0.77
0.88
0.93
1.02

Max.
Inundation
depth (m)
0.033
0.046
0.051
0.062
0.066
0.074

Table 2: Properties of debris
No.
1
2
3
4
5
6
7
8
9
10
11
12

Type

Rectangular

SphereDimensions (cm)2.5x2.5x2.5
2.5x2.5x3.0
2.5x2.5x4.0
2.5x2.5x5.0
3.0x3.0x3.0
3.0x3.0x4.0
3.0x3.0x5.0
4.5x4.5x4.5
4.5x4.5x5.0
345Mass (g)
9
10
14
17
14
18
22
47
54
21
47
75

3. Results and Discussion
3.1 Wave forces

Fig.1 Sketch of the experimental set-up

Figure 2 shows the time histories of wave forces exerted onhouse model by each cases of wave bores listed in Table 1,
without the presence of debris. Peak wave forces varied in
between 0.32 N to 11.06 N due to the wave action.

A one storey house model adopted from typical house reported
in Chen et al. (2005) was employed in this study. This model has
a dimensions of 0.15 x 0.1 x 0.05 (length x width x height in m)
and made up of wood.
2.2 Experiment cases
There were six contrasting flow conditions were generated torepresent a band of diverse inundation depths and approachingvelocities representing possible tsunami characteristics are listed

Fig.2 Wave forces on house model3.2 Debris impacts
First, time histories of wave and debris impact forces werecompared together and related maxima in each cases wereidentified. Figure 3 presents the shifts in maximum collision
forces of rectangular debris of 2.5x2.5 cross section with thelengths of 2.5, 3, 4 and 5 cm varying with the maximum
inundation depth measured near house model (at W4). Similarly,
Figures 4 and 5 show the cases of rectangular debris with thecross sections of 3.0x3.0 and 4.5x4.5 respectively. Variation ofpeak forces for spheres are exhibited in Figure 6.

Fig.3 Maximum forces induced by rectangular 2.5x2.5

As detected, rectangular debris of 2.5x2.5 cross section showsa notable increment for Cases I and II and showed quite similar
behavior of wave forces for other cases. For lengths of 4 cm and 5
cm of the same yielded a notably higher force than the wave
effects. It might be due to the increment in size and mass of theobject. Correspondingly, rectangular debris with 3.0x3.0 crosssection also behaves like rectangular 2.5x2.5 debris, where
objects with 3.0x3.0 cross section exerted comparably larger
forces than debris 2.5x2.5. However, for the 4.5x4.5 cross section
debris, the lowest flow was not able to carry the debris to thetarget location as wave force is dissipated in the object, thus force
induced on the house also reduced. However, toward the larger
flows it displayed an extensive amount of force than wave forces
as the size and mass were increased. Which implies that contact
area with the surface also important. On the other hand, spheretype of debris well behaved than rectangular debris as for eachcases additional impact forces were significantly expanded andclearly visible than the rectangular debris. It is due to their contact
area with surface is very much smaller, hence they have theability of moving even in smaller flows. In overall, additional
forces were found to be in the range of 15% to 60%. Further,
initial and impacting positions of debris has a certain role inaffecting the impact which was recognized for larger length
rectangular debris.

4. Conclusions

Fig.4 Maximum forces induced by rectangular 3.0x3.0

This study was performed to identify the impact forces andbehavior of debris brought by tsunami flow into the land area hit
a typical residential facility. Rectangular and sphere type ofdebris with different dimensions and sizes were utilized. As
observed through the experiments, debris can cause 15% to 60%
of additional impact to the structures. It was noted that, smaller
debris performed well in the smaller flows and for high flows
their impact was not significant. For larger debris, it shows aconsiderable impact forces during the higher flows and during thelower flows, sometimes it was not able to move the debris in a
desired way. Parameters such as inundation depth, distance
between debris and house, velocity, initial and impacting
positions, hydrodynamic viscosity, mass and size of debris can beconsidered as governing parameters.

References
Fig.5 Maximum forces induced by rectangular 4.5x4.5

Fig.6 Maximum forces induced by spheres

1) Chen E, Ho E, Jallad N, Lam R, Lee J, Zhou Y, Re D.D,
Berrios L, Nicolimo W and Ratti C (2005), Resettlement or
Resilience?, The Tsunami Safer House Project, International
Symposium on Disaster Reduction on Coasts ScientificSustainable Holistic, Monash University, Australia.
2) FEMA 55 (2011), Coastal Construction Manual, 4th Edition,USA.
3) FEMA P646 (2012), Guidelines for design of structures forvertical evacuation from Tsunamis, 2nd Edition, USA.
4) Murata S, Imamura F, Katoh K, Kawata Y, Takahashi S and
Takayama T (2016), Tsunami: To Survive from Tsunami, 2 nd
Edition, World Scientific Publications, Singapore.
5) Nistor I, Palermo D, Nouri Y, Murty T and Saatcioglu M(2008), Chapter 11, Handbook of Coastal and Ocean
Engineering, World Scientific Publications.

International Program of Maritime & Urban Engineering / Abstract of Master Thesis

July, 2020

Study on the effect of weld bead shape on fatigue strength of
Out-of-Plane Gusset
BUERLIHAN AYANG
Dependability and Optimum Design, Department of Civil Engineering
Key Words: Additional weld, grinding, fatigue life, out-of-plane gusset, crack initiation and propagation

1. Introduction
With the aging of social infrastructure structures, it has become
important of consideration of their maintenance. Fatigue,
corrosion are the main damage factors of steel structures. The
most important things in maintenance is the detection of damage
such as cracks and the prediction of service life. In particular,
fatigue is not limited to social infrastructure structures, but is amajor factor of damage as it is said that about 80% of damage tomachinery and structures is caused by fatigue. Therefore, it isimportant to evaluate fatigue strength correctly for properly
maintain.

2. Fatigue test

As it is noticed from Fig. 2, the additional welding withdifferent weld aim position leads to improvement of fatigue
strength of welded joints which is due to the noticeable increasingof weld toe radii and the significant decreasing of the flank angle.
Furthermore, the more fatigue strength improvements achievewith longer weld aim position for additional weld treatment
which about 1.84 times of fatigue life of AW specimen at stressrange Δσ=200MPa. On the other hand, it is also shown that thefatigue strength of B2G specimens extremely extended bygrinding method with about 6.89 times of AW specimen at stressrange Δσ=200MPa. It is needed to be mention that the different
tendency of fatigue life revealed for B0G, B2G, B4G specimens
than additional weld specimens, where no more increasing of
fatigue strength obtained for B4G specimen.

2.1 Specimens
The specimens used in this experiment were out-of-plane
gusset welded joints consisting of a same thickness for main plate
and two gusset plates as shown in Fig. 1. The parent material forthe fabrication of out-of-plane gusset was JIS G 3106 SM490YA.

Fig.2 S-N curves for fatigue life

3. Numerical analysis for fatigue tests
Fig.1 Geometry of specimens
The specimen used for fatigue tests divided to three types. As
weld (AW), additional weld (B0, B2, B4) and the additional
welding specimens treated with grinding (B0G, B2G, B4G).
2.2 Fatigue tests
Fatigue testing was performed under uniaxial tensile as thefollowing conditions: stress ratio = 0.05, loading frequency = 7
Hz, with nominal stress ranges of 125, 150, 175, 200,225 and 250
MPa. Strain gauges were affixed to 5 mm away from each boxing
weld toe on the main plate. The strain range was recorded every
5000 cycles. Crack initiation life (Nc) and Fatigue life (Nf) were
defined as the numbers of cycles at which the measured strain
decreased by 5% from the initial value and when the specimen
fractured, respectively.

3.1 FE model
3D finite element (FE) models were generated based on thespecimen geometries of the fatigue test specimens by 3D scanner.
The FE model and cross sections of weld bead for as weld andadditional weld model were illustrated in Fig. 3. Considering theshape of the test specimens and for the purpose of the reducing of
numerical costs, only 1/8 of geometry was modeled in 3D for all
types of specimen.

2.3 Fatigue test results and discussions
Fig. 2 shows the results fracture life. The horizontal axis of thegraph indicated the number of cycles, the vertical axis shows thenominal stress range, and the average design curve of the JSSC
for each fatigue class is shown by a dash line in the graph.

Fig. 3 FE model3.2 Analysis conditions and assessment method

stress intensity factor range.

The numerical analyses were conducted in two separate steps:
elastoplastic analyses and X-FEM analyses based on elastic
fracture mechanics. In addition, to evaluate the effect of the shapeof the additional weld bead on the stress concentration, the elastic
analysis was performed.
The fatigue life of crack initiation was assessed based on themethod adopted in the previous studies1). Fatigue S-S model2)
used for elaso-plastic analysis. The critical also was modified byconsidering the effect of mean stress as following equation.
𝛥𝜖𝑡𝐴(𝛔𝑚 ) = 0.305𝑁𝑐−0.580 + 0.00412𝑁𝑐−0.115
2𝐴(𝛔𝑚 ) =1
1 − 4.1 × 10−4 𝜎𝑚 − 1.6 × 10−7 𝜎𝑚2

Where σm is mean stress.
Paris’ law was adopted for the estimation of crack propagation
fatigue life. The material constants of da/dN-ΔK average design
curve were chosen based on Guidelines for Fatigue Design ofSteel Structures and Explanations by JSSC3).

Fig. 6 SIF-a relationship (Δσ = 200MPa)
In Fig. 7, the fatigue lives obtained in the numerical simulation
are compared against the experimental data for the stress range at
200 MPa.The results show that the predicted fatigue lives fairly
agree with the experimental fatigue life.

3.2 Results and discussion
The Fig. 4 shows the value of Maximum principle stress foreach type of model. The results revealed that the additional
welding treated with grinding method remarkably reduces thestress concentration.

Fig. 7 S-N curve (for AW, B4, B2G cases)

4. Conclusions
Fig. 4 The value of Maximum principle stressFig. 5 illustrates the S-S curves of elements with the highest
cumulative plastic strain for each model at stress range Δσ =200MPa. A smaller ratcheting is observed for B2G modelcompared with the AW case.

1) The local stress-strain curves, evaluated at the element with thehighest accumulation of plastic strain, show smallest sizes of thehysteresis loops and reduced material ratcheting in the models
with additional welding treated with grinding. As a consequence,
it can be concluded that this treatment is beneficial for enhancing
the fatigue crack initiation life.
2) Observing the crack propagation analyses, it can be concluded
that the additional welding treated with grinding method has alsoa small influence on crack propagation life.
3) The comparison of the numerical simulations and theexperimental data revealed that the FEA can predict with areasonable degree of agreement of the experimental fatigue life.

References1)
2)

Fig. 5 4 Local stress-strain curves (Δσ = 200MPa)
From Fig. 6, it can be found that the models with an additional
weld treated with grinding method models display overall a lower

3)

JSMS: Hirou Sekkei Binran, Yokendo Co.Ltd, (1995)
S. Tsutsumi, R. Fincato, Cyclic plasticity model for fatigue
with softening behavior below macroscopic yielding,Materials & Design, 165 (2019) 107573.
Japanese Society of Steel Construction: Guidelines for
Fatigue Design Recommendations for Steel Structures,
Gihodo Shuppan, 2012.

International Program of Maritime & Urban Engineering / Abstract of Master Thesis

July 2020

Systematic Study of Water-on-Deck Effect of Offshore Supply Vessel
Running in Stern Quartering Waves
Sreenath Subramaniam
Ship Intelligentization Subarea, Naval Architecture and Ocean Engineering
Key Words: Water-on-deck, Systematic model runs, Shallow water deck flow, Finite volume method

1. Introduction

2.2 Arrangement for measuring water level

Ships like offshore supply vessels and fishing vessels with anextended low weather deck could trap water on the weather deck
when they run in high waves. This trapped water can be critical forship transverse stability. In order to study the water-on-deck effect,a series of model experiments were conducted by Umeda, Aqmil
et al. 1) in 2015 and by Umeda, Su et al. 2) in 2016. Through theexperimental studies, a numerical model which considers thewater-on-deck effect was developed; however, the model could notpredict the actual amount of water accumulated on deck a priori,and it had to be determined by trial and error. Also, the amount ofwater accumulated on deck was not measured in the previous
experiments.
The author took up this limitation as the research aim, i.e., todevelop a numerical model which can estimate the amount of water
accumulated on deck. A new set of systematic model runs were
executed for stage-wise validation of the numerical models. New
systems were added to monitor the water trapped on deck. The
estimation of water-on-deck consists of two problems:• Water inflow and outflow through openings, based on thepressure difference.
• Water flow across the deck.
For the master's course work of the author, the second problemof water flow across the deck was solved first. A numerical model
based on shallow water flow theory was developed and validated
against the systematic model experiment data.

For measuring the varying water level in the deck, four cameras
were positioned at the four corners of the weather deck, and one
camera was placed at the forecastle with a view of the deck. Scales
were placed opposite to each camera placed on the weather deck to
detect the changing water level. In addition, four pressure sensors
were added along with the cameras. The setup is shown in Fig. 1

2. Systematic Model Runs

Cameras

Scale

SensorCameraFig. 1 Arrangement of camera and pressure sensor system
2.3 Main Observations
The author was able to confirm the presence of damping causedby the water moving across the weather deck by comparing thephase lag between ship roll motion and water motion across theweather deck, as shown in Fig. 2. The variation of damping withship speed was also confirmed through the experiment.
Fn = 0.32, 𝐻/𝜆 = 0.03, 𝜆/𝐿 =1.5
Prefilled water = 3885g

2.1 Stage-Wise Planning
In order to isolate the water-on-deck effect into distinct subphenomenon, the systematic model runs were planned in 3-stages.
Stage 1 - General Case: At first, the author executed the mostgeneral case, which is the free-running model test with all theweather deck openings, i.e., freeing ports and stern opening opened.
The results from the experiment contain the complete ship motion
and water on deck phenomenon.
Stage 2 - Evaluation of ship motion: In this set of free runs, all
the openings in the weather deck of the model ship were closed.
Thus, the results from this set of free runs contain only the ship
motion due to waves; the water-on-deck effect is not present. These
results can be used to validate the ship motion part of the code.
Stage 3 - Evaluation of ship motion and water motion on deck:
For this stage, all the openings in the weather deck were closed
similar to the previous set of model runs; however, the deck was
now filled with water in advance. Thus, the results from theexperiment contain ship motion and water motion on deck, which
can be used to validate the corresponding combined numericalcode of ship motion and deck water flow.

Phase Difference

Fig. 2 Time history of ship roll and water levels
The details of the model ship, the experimental conditions andthe observations were published in the proceedings of the 2020
spring conference of The Japan Society of Naval Architects andOcean Engineers 3).
The data from Stage 3 of the model experiments were used in thevalidation of the deck flow model.

3. Numerical simulation of deck water flow

3.2 Results of the flow model

3.1 Coordinate system and governing equations
The body-fixed coordinate system of the ship is used to describe
the motion of the fluid on the deck. When this non-inertial system
undergoes motion in the space fixed inertial system, forces are
induced on the fluid inside the deck.

The deck motion is simulated in 6 DOF using the data from thegyro, and the flow of water across the deck is numericallyevaluated using the developed model for deck water flow. The
comparison of the water heights from simulation and experimentat the location of the aft-port side sensor is shown in Fig. 4.
Fn = 0.35, 𝐻/𝜆 = 0.03, 𝜆/𝐿 =1.5
Prefilled water = 3885g
Bulwark height = 90mm

Fig. 4 Comparison of water levels between exp. and sim.
Fig. 3 Coordinate systems used to describe fluid motion
Fn = 0.35, Time = 9.18 s

The flow of water across the weather deck is physically similar
to the shallow water flow. In the governing equations, viscosity isneglected, and the forces due to ship motion appear as source terms
in the horizontal momentum equations as shown below,
𝜕(ℎ) 𝜕(ℎ𝑢)
𝜕(ℎ𝑣)
+
+
=0
𝜕𝑡
𝜕𝑥
𝜕𝑦
𝜕(ℎ𝑢)
𝜕
1
𝜕(ℎ𝑢𝑣)+(ℎ𝑢2 + 𝑔ℎ2 cos 𝜃 cos 𝜙) +
𝜕𝑡
𝜕𝑥
2
𝜕𝑦

(1)

(2)

= ℎ [−𝑔 sin 𝜃 + 𝑓𝑓𝑥 + 𝑓𝑟𝑥 ]𝜕(ℎ𝑣) 𝜕(ℎ𝑢𝑣)𝜕1++(ℎ𝑣 2 + 𝑔ℎ2 cos 𝜃 cos 𝜙)
𝜕𝑡
𝜕𝑥
𝜕𝑦
2

Deck Water

Fig. 5 Deck water flow experiment vs simulation (top-left)(3)= ℎ[𝑔 sin 𝜙 cos 𝜃 + 𝑓𝑓𝑦 + 𝑓𝑟𝑦 ]where ℎ is the free surface height, 𝑢 is the x-component of
the velocity, and 𝑣 is the y-component of the velocity. ⃗⃗⃗𝑓𝑓 is theinertial force vector arising due to the ship motion. ⃗⃗⃗𝑓𝑟 is thefrictional force between the deck and fluid, evaluated based on the
Manning’s friction slope model.4)
The above system of hyperbolic partial differential equations issolved in the discretised deck using the finite volume method. 5)
The structured mesh consists of 300 cells in the x-direction and 100
cells in the y-direction. Time steeping is done using the 3-stepstability preserving Runge-Kutta method. The intercell numericalflux is calculated using an approximate Riemann solver andreconstruction of the solution vector at cell faces from the cellaverages is carried out using forward difference together with the
Van-Albada slope limiter to avoid oscillations in the final solution.
Reflective wall boundary conditions are applied at the four wallsof the deck and time stepping is controlled using the CFL condition
to ensure stability. The Courant number is chosen as 0.9.

4. Conclusions
A new set of systematic model experiments were executed withdirect measurement of the accumulated water-on-deck. The
presence of damping was confirmed by comparing the phase lag
between ship motion and water motion 3).
A new numerical model for deck water flow based on shallow
water theory and finite volume method is developed and is
successfully validated using the systematic model experiment data.
This solves one of the problems involved in the numericalestimation of water-on-deck. The water ingress and egress problem
will be solved in the doctoral course.

References
1) Umeda, N, Aqmil, et. al.: Proc. of 15th ISSW, pp.11-16, 2016.
2) Htun, S, Umeda, et. al.: OE, Vol. 188, Article 106225, 2019.
3) Sreenath, S, Umeda, N, Maki, A, Matsuda, A: Proc. of Spring
Conference JASNAOE, pp. 623-627, 2020.
4) Tatyana, D, Alexander, K: J. Phys.: Conf. Ser. 973, 012032,
2018.
5) Fayssal, B, Mohammed, S, J of Computational and Applied
Mathematics, Vol. 234, pp. 58-72, 2010.

International Program of Maritime & Urban Engineering / Abstract of Master Thesis

July, 2020

Still Water Bending Moment Evaluation on Beam through Cargo
Loading Prediction using Kalman Filter
Rodhiatul Isnaini
Ship Structural Integrity Subarea, Naval Architecture and Ocean Engineering Department
Key Words: Kalman Filter, Cargo loading, Still Water Bending moment, Monitoring system, Inverse Problem

1. Introduction
In general, there are two types of loads that occur in a ship.
They are static and dynamic loads. A ship must be designed to beable to withstand the maximum load that is likely to be exerted on
it throughout its service time. Still water bending moment(SWBM) takes up almost 1/3 of the total bending moment and isheavily influenced by the distribution of its cargo loading, which,however, has not been investigated seriously while a lot has been
done on the wave bending moment. It may be attributed to a factthat the SWBM is regarded as a problem of “operation” rather
than “design.” It has been believed that the ship must conform to
the practice of the “good seamanship”.
Needless to say, for the more reasonable and rational ship
design, knowing the cargo loads is crucial to predict the total hull
girder vertical bending moment consisting of SWBM and wave
bending moment more accurately. Additionally, by knowing thevalue of cargo mass, local behaviour such as double bottom
bending will be able to be estimated to assure the ship safety
during operation.
Recently, the structural health monitoring (SHM) using
various sensors, is becoming more common. However, it does notmeasure the cargo loads or static loads directly and the strain
measurement data often contain noise or offset drift. Sometimes
it is improbable to utilize available sensors to get the desirable
measurements data. This research addresses these problem byusing Kalman filter and numerical model to predict the cargoloading on a beam from limited measurement data that areprovided by available types of sensor in SHM, as an inverse
problem. SWBM will be evaluated and compared to data thathave been calculated from FEM.

2. Kalman Filter Modeling1)
The ship structure can be regarded as a beam that has several
cargo holds. For example, in the case of cape size bulk carrier,there are nine cargo holds typically. The present problem can be
reduced to a prediction of the cargo mass in each cargo hold.
Kalman filter is a recursive algorithm which works based on
statistics to provides more accurate prediction of unknown
variables from a system. This algorithm combines informationfrom the system itself and a series of measurement data over time
containing noise. By combining these, the covariance of predictedunknowns is able to be minimized, therefore it yields stronger andmore accurate results from statistical perspective.

𝑥 1 𝑛+1
𝑥1 𝑛𝑤1{ ⋮ } = [𝐹] { ⋮ } + { ⋮ }𝑥 𝑁 𝑛+1
𝑥𝑁𝑛
𝑤𝑁

(1)

𝑧1 𝑛
𝑥1𝑛
𝑣1
{ ⋮ } = [𝐻] { ⋮ } + { ⋮ }
𝑧𝑀 𝑛
𝑥𝑁𝑛
𝑣𝑁

(2)

{𝑧} = {𝜀 1 , … 𝜀 𝑀 , 𝛿 1 , … 𝛿 𝑀 }(3)

Kalman filter consists of state-transition model andobservation model that are given in equations 1 and 2,
respectively. This time, cargo loading is assumed to be the statevariables, and the number of predicted state variables is denoted
as N. The observation variables, which consists of deformationand displacement data are denoted as z in equation 2 with M as
the number of sensors. It is assumed that the deformation ismeasured in terms of strain. State matrix [F] is used to link priorand posterior state variables, while observation matrix [H] is used
to correlate observation and state variables. The matrix may beconstructed by using FE model. Moreover, this research assumes
the distribution of all variables to be Gaussian with 10% white
noise represented by w and v.
Strain and displacement are measured at several points on thebeam. Strain was measured by using strain gauge. For both simply
supported and hydro-element beams (explained later), strain
gauges were placed in the same five positions. They are
represented by red dots in figure 1. Meanwhile, displacement for
simple beam was measured in nine points, while displacement ofhydro-element beam was measured in 17 points.
The displacement data are obtained by using a motion capturesystem provided by an optical tracking. This system consists ofhigh-speed tracking cameras that capture motions of an object
through specific markers (grey balls in figure 1) attached to it. The
whole process can be done by the help of motion capture software
utilization provided by the system provider.

Type

Table 1 Beam dimension
LengthWidth

Simple beam
Hydro-element beam

1000 mm1000 mm

40 mm
30 mm

Thickness
2 mm
3 mm

3. Subject Structure
3.1 Simple Beam
The first beam that is investigated in this research is simply
supported beam in air with its dimension explained further inTable 1. Experiment has been done and the output data are used
to complete prediction data for Kalman filter formulation.
This time, the beam is partitioned into twelve number of
elements, which represents the number of cargo holds. Moreover,two different number of unknown variables prediction is
performed. One being two number of unknown variable, and theother case assumes that there are nine numbers of unknown
variable. These number of predicted unknowns represents thenumber of cargo holds that are estimated. Furthermore, when two
cargo holds are assumed to be unknown, this means that 10 othercargo holds are assumed to be fixed. The same case is also applied
for the other case. This indicates that cases in hydro-element beam
analysis follow this assumption as well.
3.2 Hydro-Element Beam
The second structure is also a beam, however, floating on awater surface, (figure 1). For the modelling, a hydro-element isdeveloped. Table 1 gives further details on its dimension. This
beam is divided into 18 elements. Similar to the previous system,
two types of prediction are performed. The first one is prediction
with five unknown variables, while the second case assumed to
be predicting nine unknown variables.

Fig.1 Hydro-element beam experiment set-up

2 (a) shows the cargo loading prediction by assuming only two
unknown variables, while figure 3 (a) presents nine unknown
variables. The maximum error that is obtained from assuming two
unknown variables prediction is around 8%, while the otherprediction with nine unknowns has maximum error as big as 10%.
As for bending moment errors, they are presented in figure 2 (b)and figure 3 (b) for two unknown predictions and nine unknown
predictions, consecutively. The error that are found for bending
moment values is lower than ~10% in each node, which showsthat the Kalman filter works well in these predictions.
4.2 Prediction Results on Hydro-Element Beam
This subsection presents and discusses prediction results inhydro-element beam, including the bending moment results aswell. Identical to section 3.2, there are two outputs ofmeasurement data that were obtained from experiments. These
data are strain and displacement that were utilized in the Kalmanfilter formulation.
The first case that was being evaluated is cargo loading
prediction performed by assuming five unknowns. Figure 4 (a)and (b) shows cargo loading and vertical bending momentcompared to its reference, consecutively. The error that have been
calculated from this case are lower than 10% for both cargoloading prediction and bending moment, depicting that the filter
works well for this case.

4. Results and Discussion
4.1 Kalman Filter Implementation on Simple Beam(a)(b)
Fig. 2 (a) Two cargo loading prediction (blue) in comparison withthe reference value (red) (b) Vertical bending moment
obtained from prediction (blue) in comparison with its
reference value (red)(a)(b)
Fig. 4 (a) Five cargo loading prediction (blue) in comparison withthe reference value (red) (b) Vertical bending moment
obtained from prediction (blue) in comparison with its
reference value (red)The second prediction is carried out by assuming there arenine unknown variables. This arrangement, however, resulted instate variables prediction with big errors compared to its reference.
This leads to high error in still water bending moment results aswell. The maximum error for the bending moment result reach
50%, which is not an ideal value for error.

5. Conclusions
In this research, Kalman filter is utilized as a tool to develop
a monitoring system. The system has an ability to predict cargoloading with different cases of unknown variables, as long as
measurement data provides values with acceptable error.
(a)(b)
Fig. 3 (a) Nine cargo loading prediction (blue) in comparison withthe reference value (red) (b) Vertical bending moment
obtained from prediction (blue) in comparison with its
reference value (red)Prediction using Kalman filter based on strain anddisplacement data from experiment has been investigated. Figure

References
1)

Perisic, N. and Kirkegaard, PH.: Low-Cost Tower Root
Fatigue Load Estimation for Structural Health Monitoring of
Grouted Connections in Offshore Wind Turbines, Proc. KeyEngineering Materials, 2013

International Program of Maritime & Urban Engineering / Abstract of Master Thesis

July, 2020

Development of Iso-Geometric Analysis Software JWRIAN-IGA
for Weld Bead Geometry Design Using Bayesian Optimization
Thein Lin Aung
Department of Naval Architecture and Ocean Engineering
Key Words: Isogeometric Analysis, Geometry optimization, Fillet weld bead, Bayesian

1. Introduction

3. Methodology

Iso-Geometric Analysis (IGA) is an emerging method ofnumerical simulation which uses the same data of Computer Aided
Design (CAD) for both geometry description and structure analysis.
In CAD geometry, B-Splines or Non-Uniform Rational B-Splines(NURBS) are commonly used to describe free-form surfaces suchas the body of an automobile or the hull of a ship. On the other
hand, traditional Finite Element Analysis (FEA) typically uses
Lagrange polynomials to represent geometry and solution fields.
Because of this difference, pre-processing such as meshing ofgeometry models take up much of engineering analysis time1).
Moreover, it is difficult to feed back the analysis results to CAD
environment for further design modifications. IGA is a unified
method which hopes to remedy these problems by integrating CAD
and FEA. This is achieved by using the same basis functions torepresent the geometry and solution fields. In this way, the analysis,modelling (such as meshing by pre-processing) and evaluation(such as visualization by post-processing) become seamless andefficient, because they are working on the same platform.
In shape optimization, the parametrization of the model for thenumerical evaluation of the objective function and exposing the
optimization variables are of paramount importance2). NURBS
geometries provide a convenient way to vary shape parameters
since control points can directly be used and high order continuity
ensures the smoothness of the boundary. Moreover, IGA enables
the use of modified geometry for simulation without further preprocessing which makes the process highly automatic. In addition,
the optimized geometry can be exported to CAD since the final
output is also a NURBS geometry.
Bayesian optimization is a stochastic optimization method basedon Bayes’ Theorem in probability theory. It determines the
sampling parameters by approximating the objective function fromobserved data by Gaussian Process. It strives to find the optimumin the fewest number of objective function evaluations as possible.
Hence it is suited for applications where the function evaluation iscomputationally expensive such as FEA simulations.
In this study, IGA software JWRIAN-IGA is developed and theshape of a fillet weld bead under transverse load is optimized toreduce stress concentration with minimum weld volume using
Bayesian Optimization.

In Isogeometric Analysis, B-Splines or Non-Uniform Rational
B-Splines (NURBS) functions are used as basis functions (shapefunctions) to represent both geometry and solution fields. B-Spline
basis functions of order 𝑝 are defined by Cox-de Boor recursion
formula3):
𝑁𝑖,0 (𝜉) = {𝑁𝑖,𝑝 (𝜉) =1 , if 𝜉𝑖 ≤ 𝜉 < 𝜉𝑖+1
and
0 , otherwise(1)
𝜉𝑖+𝑝+1 − 𝜉
𝜉 − 𝜉𝑖𝑁𝑖,𝑝−1 (𝜉) +𝑁𝑖+1,𝑝−1 (𝜉)𝜉𝑖+𝑝 − 𝜉𝑖𝜉𝑖+𝑝+1 − 𝜉𝑖+1

where ξ is the parameter value and 𝜉𝑖 denotes 𝑖𝑡ℎ knot value
in the knot vector {𝜉1 , 𝜉2 , … , 𝜉(𝑛+𝑝+1) }. A knot vector is a nondecreasing set of real numbers in the parameter space which define
the mesh boundaries.
NURBS basis functions are calculated by adding an additional
weight parameter 𝑤𝑖 to each control point and taking theweighted average of B-Spline functions as shown in Eq. (2).
𝑝

𝑅𝑖 (𝜉) =𝑁𝑖,𝑝 (𝜉)𝑤𝑖𝑛∑𝑖=1 𝑁𝑖,𝑝 (𝜉) 𝑤𝑖

The symmetry model of a fillet welded T-joint is shown in Fig.
1. The geometry of fillet weld along with the base and web plates
were created as separate NURBS solid patches. Second orderfunctions were used for all patches. Each control point has 3translational DOFs and no rotational DOF. Gaussian quadrature
with full integration points was used for the integration in thedomain. Multi-patch coupling is done using shared control points
at patch boundaries. Linear constitutive relation is assumed with
𝐸 = 210 GPa and 𝑣 = 0.3 for the whole model.

2. Objectives
The objectives of this thesis are•To find the optimal shape of fillet weld using solid multipatch IGA and Bayesian Optimization•
To study the feasibility and accuracy of IGA and Bayesian
Optimization in geometry optimization problems•
To identify improvements and future research opportunities
in IGA shape optimization.

(2)

Fig. 1 Base geometry of the fillet weld model with patch
numbers

The positions of control points on the top surface of fillet weld
were used as variables to control the shape of the weld bead asshown in Fig. 2. After the control points were modified using theparameters from optimizer, the geometries are then refined usingh-refinement (knot insertion) process to get control point mesh andknot-span mesh required for the numerical simulation. The
geometry was not affected by the refinement process. The objectivefunction was defined as a combination of maximum equivalent(Von-Mises) stress and weld volume as shown in Eq. 3.
min: 𝑓(𝒅) =𝑚𝑎𝑥⁡(𝜎𝑒𝑞 )
𝜎𝑒𝑞𝑜

+ 𝜆.

𝑉
𝑉𝑜

(3)
Fig. 3 Reduction of cost in each iteration step.

The normalization constants were set to 𝜎𝑒𝑞𝑜 = 100 MPa and𝑉𝑜 = 7000 mm3, and the penalization parameter was set to 𝜆 = 1.
The optimization variables were bounded by their respective
extrema and also subject to linear constraint relations to prevent
infeasible geometry configurations. The weld shape was assumed
to be uniform in the longitudinal direction. The total number of
optimization variables was 10 and number of iterations was 30.
Isogeometric Analysis was used as a black-box function tocompute equivalent stress and weld bead volume in each geometry
iteration. Newly developed in-house code JWRIAN-IGA was used
for IGA and open source code BoTorch4) was used for Bayesian
Optimization.

Fig. 2 Parametrization of weld geometry with optimization
DOFs shown in red arrows.

4.

Results and Discussion

The values of the objective function (cost) for each iteration areplotted in Fig. 3. The cost decreased significantly during the firstfew iterations where random parameter sets were used to fit theGaussian model. After that, the cost generally decreased but
sometimes increased in a few iterations because of the exploratory
part of global optimization. In those cases, only current best valueof the cost was plotted. Overall, it can be concluded that the
optimization algorithm was effective since it required only 30
iterations to optimize for 10 parameters.
The shapes of weld bead for some iterations are shown in Fig. 4.
The surface of the weld was a second order B-Spline curve andhence it was always smooth regardless of the discretization bycontrol points. The numerical analysis also benefited from theexact representation of weld geometry because the shape remained
exactly the same before and after mesh refinement. Generally, boththe maximum equivalent stress and weld volume decreased as theiteration number increased. Each function evaluation tookapproximately 5 minutes and the total run time was under 3 hours
by sequential iterations. This duration can be made shorter by
employing parallel simulations during optimization.

Fig. 4 Shapes of weld bead in selected iterations.

5. Conclusions
Isogeometric Analysis and Bayesian Optimization were used todesign the optimal shape of a fillet weld to minimize stressconcentration and volume of weld metal. From this study, thefollowing conclusions can be drawn.
1. Using NURBS control points for weld shape resulted insmooth geometry parametrization
2. IGA method allowed the use of exact geometry for both
optimization and analysis models
3. Bayesian Optimization achieved optimal weld shape in only
30 iterations with total run time of less than 3 hours

References
1)

2)
3)
4)

Hughes, T. J. R., Cottrell, J. A. & Bazilevs, Y. Isogeometricanalysis: Comput. Methods Appl. Mech. Eng. 194, 4135–
4195 (2005).
Frazier, P. I. A Tutorial on Bayesian Optimization. (2018).
Piegl, L. & Tiller, W. The NURBS Book. (Springer Berlin
Heidelberg, 1995). doi:10.1007/978-3-642-97385-7.
Balandat, M. et al. BoTorch: Programmable Bayesian
Optimization in PyTorch. (2019)

International Program of Maritime & Urban Engineering / Abstract of Master Thesis

July 20, 2020

Study of Automatic Berthing by Evolution Strategy with CovarianceMatrix Adaptation (CMA-ES) in Real Time
Dimas M. Rachman
Ship Intelligentization Subarea, Naval Architecture and Ocean Engineering Department
Key Words: Automatic berthing, CMA-ES, Path-following, Optimal control theory

1. Introduction

𝛿(𝑡)2 − (𝛿𝑚𝑎𝑥 )2 ≤ 0, 𝑛(𝑡)2 − (𝑛𝑚𝑎𝑥 )2 ≤ 0Over the past few decades, researchers have applied various
mathematical strategies in designing a berthing controller. As early
as 1987, Koyama et al.1) pioneered the first systematic study about
automatic berthing in which it is modeled as a linear quadratic
regulator (LQR) problem where the collision risk with berth’s wallis also considered. Shouji et al.2) attempted to use nonlinear,constrained optimal control theory in the berthing controller. Theymodeled the berthing operation as a minimum-time problem whichmeans, the ship shall arrive at the berthing point within the required
minimum time. In solving the optimal control problem, theyutilized a numerical technique called sequential conjugate
gradient-restoration algorithm (SCGR).
More about optimal control problem, Maki et. al.3) are among
the first to investigate the use of CMA-ES for automatic berthing
where they utilized CMA-ES to generate optimal trajectories ofberthing operation along with the optimal control. They modeled
the automatic berthing operation as a minimum-time problem withconsideration of collision risk to the berth’s wall. The study showsthat CMA-ES successfully solves the nonlinear berthing problem
despite the lengthy computational time. To further investigate its
use in the field of naval architecture, the author attempts to employ
CMA-ES in the making of an automatic berthing controller thatworks in real time. Since CMA-ES is a lengthy calculation, it
seems not suitable to solve a difficult optimal control problem inreal time like automatic berthing. However, with a few
considerations in the problem definition, it is plausible to use
CMA-ES in real time, which is the main scope of this thesis.
Provided that 3) had established an effective means to generate
optimal trajectories, the author adopts a path-following approach,where the controller will try to control the ship to follow thereference trajectories (i.e. optimal trajectories). Further, coursechanging experiment is conducted to see CMA-ES performance.

𝑛(𝑡) = |𝑛(𝑡)| if 𝑡 ≤ 𝑡𝑠

2. Generation of Reference Trajectories
The optimal control problem for generating reference
trajectories can be stated as follows,
Determine𝛿 ∗ , 𝑛∗ , 𝑡𝑓 , 𝑡𝑠

thatminimize

𝐽𝑟𝑒𝑓 = max (‖𝐱(𝑡𝑓 ) − 𝐱𝑓𝑖𝑛 (𝑡𝑓 )‖ , 𝐱𝑡𝑜𝑙 (𝑡𝑓 ))𝑛(𝑡) = −|𝑛(𝑡)| if 𝑡 > 𝑡𝑠 ,where 𝐱 = [𝑥0 𝑢 𝑦0 𝑣𝑚 𝜓 𝑟]𝑇 is the state vector, and𝐮 = [𝛿 𝑛]𝑇 is the control vector. 𝛿 ∗ and 𝑛∗ are the optimal
rudder angle and optimal propeller revolution number, respectively.

3. CMA-ES
Evolution strategy with covariance matrix adaptation (CMAES) is a state-of-the-art optimization technique that can guarantee
global optimum at the cost of computational time. In solving the
optimal control problem, CMA-ES will generate a population of
the solution and calculate its mean along with the covariancematrix. After that, the samples of population are ranked based on
the objective function. Those which are not feasible are eliminated.
From the remaining samples, the evolution paths (direction of thenext generation of population) is calculated. At the same time,covariance matrix is updated, however, with the mean from thepreviously generated population (both feasible and infeasible
population). This covariance matrix is used to determine the step
size for the next generation. Finally, a new population is generated
based on the updated covariance matrix, the updated step size, andthe evolution path. This process will continue until a convergencecriterion is met.

4. Real Time Simulation of Automatic Berthing
CMA-ES will be used to solve a path-following optimal control
problem. A path-following is a problem where the state vector attime 𝑡, or 𝐱(𝑡), shall be very close to the reference state vector at
the same time, or 𝐱𝑟𝑒𝑓 (𝑡) . However, it is more convenient
numerically to discretize the problem into multiple simpler
terminal-control problems, i.e., solve for the optimal control at atime 𝐮∗ (𝑡) that minimizes the difference between the state vector
at a future time 𝐱(𝑡 + Δ𝑡) and the reference state vector at a future
time 𝐱𝑟𝑒𝑓 (𝑡 + Δ𝑡).

2

𝑡𝑓

2

∫ ‖𝐱(𝑡) − 𝐱𝑓𝑖𝑛 (𝑡𝑓 )‖ 𝑑𝑡
04𝑡𝑓+ 𝑤 ∑ ∫ |min[(𝑦𝑐,𝑖 − 𝑦𝑤𝑎𝑙𝑙 ), 0]| 𝑑𝑡
𝑖=1 0

subject to𝐱̇ (𝑡) = 𝐟(𝐱(𝑡), 𝐮(𝑡), 𝑡)𝐱(0), 𝐱𝑓𝑖𝑛 (𝑡𝑓 )Fig. 1 Path-following terminal control problem

5. Course-Changing Experiment with CMA-ES

The objective function for this path-following problem is2

𝐽𝑝𝑓 = ‖𝐱(𝑡 + Δ𝑡) − 𝐱𝑟𝑒𝑓 (𝑡 + Δ𝑡)‖ .

(1)

The prediction time Δ𝑡 is taken as 0.8s, as this is the time that isrequired for the rudder of the Esso Osaka scale model ship to rotatefrom the maximum negative angle to the maximum positive angle.
There are six reference trajectories to be followed. It is can beseen that CMA-ES successfully controlled the ship to follow thereference trajectories under constant wind disturbance. An
example is shown in Fig. 2 below. Bang-bang control (frequentswitch from maximum to minimum control) happens almost for allreference trajectories. This is because we divide the problem into
multiple terminal control problems which are independent to eachother, so minimizing bang-bang is one of the works in progress.

As a first step to put CMA-ES on practice, it is utilized to control
the Esso Osaka scale mode ship in a course-changing experiment.
The objective function for this experiment is similar to that of aone-state path-following problem. It can be expressed as
𝐽𝑐𝑐 = (𝜓(𝑡 + Δ𝑡) − 𝜓𝑓𝑖𝑛 )2 .

(3)

The experiment was done with a constant propeller revolution
number, with three desired courses 𝜓𝑓𝑖𝑛 throughout theexperiment (-20 from 0 s to 60 s, +20 from 60 s to 95 s, and 0 until
the end). Despite the bang-bang control of rudder angle, CMA-ES
proved itself worthy to control the ship in real time and in realsituation at the Inukai pond without the need to execute parametertuning like in classic PID controller.

Fig. 4 Rudder angle and yaw angle in the experimentFig. 2 Example of path-following simulation with CMA-ES,
left: in red, the controlled ship; right: the controls.
However, sometimes when the trajectories are difficult to follow,
CMA-ES with the objective function (1) failed to control the ship,especially for yaw and its rate. One of the ways to overcome this,is by modifying the objective function with a weight matrix 𝐖 =[1 1 1 1 10 10]𝑇 to cover the uncertainties and to give asense of importance to yaw and its rate, since they are difficult tocontrol. The weighted objective function is expressed as
𝐽𝑝𝑤 = [𝐱(𝑡 + Δ𝑡) − 𝐱𝑟𝑒𝑓 (𝑡 + Δ𝑡)]𝑇

𝐖 [𝐱(𝑡 + Δ𝑡) − 𝐱𝑟𝑒𝑓 (𝑡 + Δ𝑡)].

(2)

Fig. 3 below shows CMA-ES controlled the ship to follow thereference trajectory by using the objective function shown in (2).

6. Conclusion
CMA-ES is known as the state-of-the-art optimization technique
but with long computational time. By simplifying the optimal
control problem, the computational time can be reduced so that it
can be used in real time (almost instantaneous). This has been
investigated in this thesis, under constant wind disturbance,
numerically. Further, as an initial step to apply CMA-ES in reality,
a simple course-changing experiment was done. In this experiment,
CMA-ES successfully controlled the ship to the desired courses
and gave a good performance under real random external
disturbances at the Inukai pond.
An important aspect is to avoid using weight matrix to create amore robust controller, i.e., no need to modify the objectivefunction to comply with certain trajectories. One of the ways is bynormalization of the state variables, which is a work on progress.

References

Fig. 3 Simulation with weighted objective function, upper: in
red, the controlled ship; bottom: the controls.

1) T. Koyama, J. Yan, and J. K. Huan, “A systematic study on
automatic berthing control (1st report),” Journal of the Societyof Naval Architects of Japan, vol. 1987, no. 162, pp. 201–210,
1987.
2) K. Shouji, K. Ohtsu, and S. Mizoguchi, “An automatic berthing
study by optimal control techniques,” IFAC Proceedings
Volumes, vol. 25, no. 3, pp. 185–194, 1992.
3) A. Maki, N. Sakamoto, Y. Akimoto, H. Nishikawa, and N.
Umeda, “Application of optimal control theory based on theevolution strategy (CMA-ES) to automatic berthing,” Journalof Marine Science and Technology, vol. 25, no. 1, pp. 221–233,
2020.