Development of a Semi-analytical Dynamic Thrust Force Model


 A moving water mass generates force which is exerted on its moving path. Cyclone generated storm surge or earthquake generated tsunami are specific examples of moving water mass the generates force along the coasts. In addition to human lives, these moving water masses cause severe damages to the coastal infrastructure due to tremendous force exerted on these structures. To assess the damage on these infrastructures, an essential parameter is the resultant force exerted on these structures. To evaluate the damages, there is hardly any quantitative method available to compute this force. In this paper we have developed a semi-analytical model, named as Dynamic Force Model (DFM), by using Variational Iteration Method to compute this force. As governing equations, we have used the Saint Venant equations which are basically 1D shallow water equations derived from the Navier-Stokes equations. The verified, calibrated and validated DFM is applied in Bangladesh coastal zone to compute dynamic thrust force due to tropical cyclone SIDR.


INTRODUCTION 26
A moving water mass generates force due to momentum created by the water mass itself 27 and the acceleration due its movement. Specific examples of moving water mass are surge 28 waves generated by cyclones and tsunami waves generated by earthquakes that strike the 29 coast with considerable magnitude of force. Though the sources of these moving water 30 masses (for example storm surge and tsunami) are different, the physical characteristics of 31 their wave propagation (in deep water), nonlinear transformation (in shallow water) and 32 runup (land inundation) are identical (Yim, 2005). The rapid changes in hydrodynamics in 33 near shore substantially damage infrastructures of coastal regions. So, from the major 34 concepts of coastal engineering, storm loading (Clark, 1982) (wind load, waves, hydrostatic 35 and hydrodynamic loads) and tsunami induced forces are necessary in a coastal building 36 code. Water related loads accompanying a coastal storm are often not adequately addressed 37 in coastal building code, with wind generated/ induced waves producing the most critical 38 and complex forces to which the coast and its structures are subjected (Clark, 1982). Few 39 current building codes of coastal structures such as the city and county of Honolulu Building 40 Code (CCH); the 1997 Uniform Building Code (UBC 97); the 2000 International Building 41 Code (IBC 2000), the SEI/ASCE 7-02 (ASCE 7), and the Federal Emergency Management 42 Agency Coastal Construction Manual (FEMA CCM) incorporates the effects of 43 hydrodynamic loads (surge/tsunami induced) (Pacheco and Robertson, 2005). The 44 generalized equations of contributing force components of these codes are discussed in Yeh 45 and Palermo, 2009;Yim, 2005. In these codes, computation of hydrostatic 46 force is the multiplicative combination of gravitational acceleration and square of specific 47 energy. Hydrodynamic force is calculated as the drag force acting on the direction of the 48 uniform flow. But these codes have over-simplified the process of estimating the 49 hydrodynamic loads (Yim, 2005). The forces are calculated for an instantaneous velocity of 50 moving water mass for design purposes of coastal infrastructures. But these codes neglect 51 the non-linearity of the variation of the total force generated during a cyclone or Tsunami. reference to a number of cases where signs of under-seepage piping erosion is thought to be 68 responsible for the failure of the levees (Seed et al. 2006). Wolff (1997) analyzed that the 69 probability of levee failure due to under-seepage is among the main causes of levee failures 70 and increases as the floodwater elevation increases. Piping may occur because of the 71 presence of seepage flow through weak layers, cracks, structural joints, tree roots, and 72 animal burrows in the embankment, without the necessity for the water level to have reached 73 the full height of the embankment (Foster et al. 2000;Fell et al. 2003; Richards and Reddy 74 2007). Rapid rise in water level and the prolonged immersion in high water also facilitate 75 embankment failure during storm surge (Luo et al, 2016). The former induces a tremendous 76 impact effect, which significantly decreases the resisting forces against seepage failure by 77 increasing erosion power of seepage flow (Luo et al., 2013(Luo et al., , 2016. Extensive research has 78 been carried out focusing on seepage failure mainly focuses on steady flow condition 79 (Sellmeijer, 1988, Sellmeijer and Koenders , 1991, Asaoka and Kodaka, 1992 2008, Maknoon and Mahdi, 2010) and transient flow condition (Ozkan, 2003 (Adomian, 1983(Adomian, , 1988) and Liao's homotopy analysis method 106 (Liao,1997). With these methods, most PDEs can be approximately solved without 107 linerization or weak linerization or small perturbations. However, the approximation 108 obtained by Adomian's method could not always satisfy all its boundary conditions, leading 109 to error near boundaries. A successful approximation of solution for partial differential 110 equations is established with no boundary problems by Ji-Huan He (He, 1997(He, , 1998(He, , 1999a 1999b) which is known as Variational Iteration Method (VIM). This method is based on 112 analytically solving the equations that result from discretizing the spatial coordinates of a 113 partial differential equation (Vidts and White, 1992). 114 Governing equations used in this study are the Saint-Venant equations derived from Navier-115 Stokes equations with one dimensional shallow water approximation. The equations are 116 characteristically non-linear partial differential equations. Functional relations for some of 117 the variables in the governing equations are not available. So, the equations are solved semi-118 analytically by applying coupled VIM and discretization method. In this study, this coupled 119 solution method is termed as 'semi-analytical'. 120

MODEL DEVELOPMENT 121
Assumptions to the development of Saint-Venant Equations 122 In Fluid Dynamics, the Saint-Venant equations (Raisinghania, 2003) were formulated in the 123 19th century by two mathematicians, Adhémar Jean Claude Barréde Saint-Venant and 124 Bousinnesque (Raisinghania, 2003). Saint-Venant equations are derived from Navier-125 Stokes equations (Vreukdenhil, 1994) for shallow water conditions (Kubatko, 2005).The 126 Navier-Stokes equations represent a general model used to model fluid flows (Bessona et 127 al., 2007;Bulatov, 2013). On the other hand, a general flood wave for 1-D situation 128 (Bessona et al., 2007;Bulatov, 2013) can be described by the Saint-Venant equations. 129 During development of the Saint Venant equations, following assumptions (Vreukdenhil,130 1994) were made: 131 • Flow is unsteady and non-uniform. 133 • Hydrostatic pressure prevails and vertical accelerations are negligible. 134 • Streamline curvature is small. 135 • Bottom slope of the channel is small. 136 • The fluid is incompressible. 137 • The gravity force is the only force which is taken into account. So, the influence 138 of the Coriolis force is neglected. 139 140 In addition to the above assumption, following assumptions are made for this study: 141  The steady uniform state of flow is expressed by Manning's equation which is used 142 as the initial condition for solution. 143  Manning's coefficient n is used to describe the resistance effects. 144 These two additional assumptions do not affect the basic governing equations. 145

Conservation of mass 146
Law of mass conservation says that mass can neither be created nor destroyed. When this 147 law is applied in any control volume of fluid (water), it means net rate of change of mass is 148 equal to the net change of mass due to inflow and outflow. 149 For fluid dynamics, the continuity equation (Vreukdenhil, 1994) becomes 150 Where h is the water depth and u is the flow velocity. 152

Conservation of momentum 153
This law states that the rate of change of momentum (mass times acceleration) in the control 154 volume is equal to the net forces acting on the control volume. 155 Inside a control volume, conservation of momentum can be expressed (Vreukdenhil, 1994) 156 as: 157 The pressure force is thus equal to: 173 and is a function of how filled the channel is. In other words, it is a function of depth h. 175 Taking the h derivative (which will be needed later), we have: 176 The gravitational force is the weight of the water slice projected along the x-direction, which 179 is mg (= ) times the sine of the slope angle , 180 The frictional force is the bottom stress multiplied by the wetted surface area 182 The bottom stress is proportional to the square of the velocity. Invoking a drag coefficient 184 Cd, we write: 185 The gradient of the pressure force becomes 196 In this equation the ratio of the cross-sectional area A over the wetted perimeter p, which 200 has the dimension of a length, is defined as 201 This is called hydraulic radius. Because we consider wide water body, which is much wider 203 than the depth, the wetted perimeter is generally not much more than the width ( ≈ ), 204 and so the hydraulic radius is approximately 205 In present study, all these equations are solved semi-analytically to compute spatial variation 216 of forces exerted by the wind and water mass. 217 Variational Iteration Method (VIM) 218 In 1978, Inokuti (Inokuti, 1978) proposed a general Lagrange multiplier method to solve 219 non-linear problems, which was first proposed to solve problems in quantum mechanics 220 (Inokuti, 1978). The main feature of the method is as follows: the solution of a mathematical 221 problem with linearization assumption is used as initial approximation or trial-function, then 222 a precise approximation at some special point can be obtained. Considering the following 223 general non-linear system, 224 where L is a linear operator, N is a non-linear operator and g(t) is a known analytical 226 function. 227 Above method was modified into an iteration method (He, 1997;Finlayson, 1972;He, 228 1998;He, 1999a; Mungkasi and Wiryanto, 2015) in the following way: 229 where is a general Lagrange multiplier (Inokuti, 1978), which can be identified optimally 231 via theVariational theory (Inokuti, 1978;He, 1999b;Nayfeh, 1985;Finlayson, 1972), the 232 subscript n denotes the n th approximation, and un is considered as a restricted variation 233 (Finlayson, 1972), i.e. ̅̅̅ = 0 234 The method is shown to solve effectively, easily, and accurately a large class of non-linear 235 problems with approximations converging rapidly to accurate solutions. 236 VIM is now widely used by many researchers to solve linear and nonlinear partial 237 differential equations. The method introduces a reliable and efficient process for a wide 238 variety of scientific and engineering applications, linear or nonlinear, homogeneous or non-239 homogeneous, equations and systems of equations as well. It was shown by many authors 240 (He, 1997(He, ,1998(He, , 1999a(He, , 1999bAbdou and Soliman, 2005) that this method is more 241 powerful than existing techniques such as the Adomian method (Adomian, 1988;Hagedorn, 242 1981) or perturbation method (Holmes,2009 Moreover, the power of the method gives it a wider applicability in handling large number 250 of analytical and numerical applications. 251

Solution of the governing equation 252
Equations (14) and (17) which are known as Saint Venant equations can be re-written as 253 And 255 Where is air density, is water density, is cyclonic wind speed, u is water velocity, 257 Cd is water drag coefficient and Cw is wind drag coefficient. 258 We express the steady uniform state of flow by Manning's flow equation (Chow, 1959) as 259 an initial condition for the Equations (20) and (21), which means 260 where n is Manning's roughness coefficient and Rh is hydraulic radius of the channel. 262 The correction functional for Equation (21) is 263 The stationary condition is given by 265 Equation (24) The zeroth approximations   is the steady and uniform condition. 279 Using Equation (29), we get the 2 nd iterated formula as 281 Using definite integral, we have 291 Equation (34) If we do not consider any wind term in the momentum equation, the solution becomes 301  Let u(x,t) denotes the water velocity (computed by using water depth and wind speed as 307 input) at spatial location x and temporal state t. 308 If a water particle travels along the curve x(t) through the water body, the rate of change of 309 velocity creates two kinds of accelerations: (a) Local acceleration and (b) Convective 310 acceleration 311 The acceleration for water particle can be calculated by using where, is called local acceleration i.e temporal variation and is called convective 314 acceleration when the particles move through regions with spatially varying velocity. 315 Since the water under consideration is moving, it is acted upon by external forces which will 316 follow Newton's second law. i.e. ma F  317 Therefore, the force for moving fluid can be calculated as 318 which can also be considered as the dynamic thrust force due to cyclone generated surge 320 wave. 321 Using Equation (35) Again, using Equation (35), we can compute the convective acceleration as

(40) 327
Now we can write Equation (38) as Equations (41) To verify the DFM algorithm and code, the velocity field computed by Equation (35)  In Case-2, water depth is varying along the channel, but no wind is blowing (Figure-4) The transfer of momentum from air to the ocean surface produces a wave field that imposes 463 shear stress which ultimately causes the rise of water level (Powell et al, 2003). The drag 464 coefficient Cw controls the transfer of momentum from air to water surface in the governing 465 momentum equation (Donelan, et al., 1996). It mainly depends on wind speed. In modeling 466 storm surge, this coefficient is a crucial parameter (Drews, 2013). Past studies show that 467 this co-efficient either remains constant or increases with the increase of wind speed 468 (Charnock, 1955;Large and Pond, 1981;Smith and Banke,1975;Bender et al.,2007;Kim 469 et al.,2008;Kohno and Higaki, 2006;Zhizhua et al., 2010). In the present study, it is found 470 that wind drag co-efficient has a non-linear relation with the wind speed (see Figure-8 and  471 Equation (50)). 472

MODEL VALIDATION 473
Model validation means ability of the model to represent the real world for which the model 474 is developed (Trucano et al., 2006). As real-world data of thrust force for any moving water 475 mass is not available, validation of DFM is performed from two different perspectives. 476 These is several times higher than storm surge (water depth above sea-bed can be as high as 4000m, 497 whereas, tidal range of a storm tide can vary between 6m to 7m depending on the tidal 498 amplitude and strength of cyclone). After reaching the coast, both the tsunami wave and 499 storm surge strikes the coastal structure with a force which is directly proportional to the 500 energy (Venkoba, 2017). This means that after reaching the coast, the water mass created 501 due to tsunami and moving with a certain velocity will exert more force on a coastal 502 structure compared to a moving water mass moving with the same velocity but created due 503 to cyclonic storm surge. The reason isthe moving water mass created by tsunami contains 504 more energy compared to the moving water mass created by storm surge (Dutykh  For storm surge, the velocity-force curve is generated by applying DFM in the hypothetical 523 channel which was used for model verification (Figure-1). In this case, we used case-3 524 model setup (Figure-6), but to get a reasonable range of force to generate velocity-force 525 curve, we have used a higher range of wind speed compared to the wind speed range shown 526 in Figure Comparison of velocity-force curve between tsunami and storm surge is shown in Figure-531 9. The results show that for the same velocity, tsunami force is several times higher than 532 storm surge force. With increasing velocity, this difference also increase. We have already 533 discussed why tsunami force should be higher than storm surge force. The velocity-force 534 curve in Figure-9 shows that DFM realistically computes force due to cyclonic storm surge. To quantify the visual qualitative color comparison, a one-to-one function is defined among 551 the colors between the two maps. In this way, it is possible to quantify the overlapped color 552 and deviated color between the two maps. The overlapped color shows the similarity, and 553 the deviated color shows the dissimilarity between the maps. The comparison is shown in 554 The DFM model has numerous prospects in improving coastal resilience in particular and 562 structural resilience in general. The model is not only able to compute dynamic thrust force 563 due to cyclone generated storm surge, but also can compute the dynamic thrust force exerted 564 by any moving water mass on a structure in its flow path (for example due to tidal wave, 565 wind wave, flows in ocean, estuaries and rivers, tsunami etc.). The output from the model 566 can be used as one of the design criteria of infrastructures on which force is exerted by the 567 moving water mass. As the DFM is capable of calculating spatial and temporal changes of 568 thrust force, it can also be considered as one of the monitoring parameters of extreme 569 climatic conditions related to moving water. As an illustration, the model is applied to 570 compute dynamic thrust force due to tropical cyclone SIDR. 571 572 Tropical 2016). Simulated thrust force due to cyclone SIDR at the time of landfall is shown in Figure-577 12. Table-2 shows the maximum thrust force at different coastal districts in Bangladesh (the 578 district names are shown in Figure-12) due to cyclone SIDR. 579 Bhola (see Figure-12) are impacted with very high magnitude of thrust force (magnitude 584 greater than 50 kN/m, see Table-2). In the post-disaster damage and loss assessment report, 585 these districts are also identified as the worst affected coastal districts in terms of economic 586 and infrastructure loss (GoB, 2008). Due to anti-clockwise rotational impacts of cyclones in 587 northern hemisphere, locations on the right side of cyclone tracks are the most affected. 588 According to DFM simulation, coastal districts of Patuakhali, Barguna and Bhola are the 589 worst affected regions which are situated at right side of the cyclone track and at the same 590 time close to the path of cyclone track. On the other hand, although Pirojpur, Jhalokati and 591 Bagerhat are situated close to the cyclone track (see Figure-12), magnitude of thrust forces 592 are relatively less in these districts (less than 50 kN/m) because these districts are located at 593 the left side of the cyclone track. The red lines in Figure-12 shows the coastal embankments 594 (locally known as polders) to protect the lands from surge wave. DFM simulated thrust force 595 inside the protected land is only due to cyclone wind, not due to combined impact of cyclone 596 wind and moving surge. But the thrust force in the unprotected land is due to combined 597 impacts of moving surge and cyclone wind. Thrust force simulated by DFM can be applied 598 to assess probability of failure of embankments during propagation of surge wave and 599 possible damage of households both inside and outside the protected land. These data can 600 also be used to further strengthen these structures or during construction of new structures 601 in the coastal region. 602 increasing non-linearity in the system, which is an inherent drawback of numerical 618 solutions. This shows model algorithm and code are appropriately implemented. To apply 619 the model in simulating dynamic thrust force due to cyclone generated storm surge, 620 Manning's n and wind drag coefficient are used as the calibration parameters. For the DFM, 621 an empirical equation is developed by relating wind drag co-efficient with wind speed which 622 can be used for any specific application (related to wind driven water flow) of the model. It 623 is found that wind drag coefficient increases non-linearly with the increase of wind speed. The views expressed in this work are those of the creators and do not necessarily represent 637 those of DFID and IDRC or its Board of Governors. Website: www.deccma.com. 638

ETHICAL STATEMENT 639
The manuscript submitted in the Journal of Natural Hazard, Springer is entitled 640 "Development of A Semi-Analytical Dynamic Thrust Force Model". It has not been 641 published elsewhere and that it has not been submitted simultaneously for publication 642 elsewhere. 643