Higher-order fem for nonlinear hydroelastic analysis of a floating elastic strip in shallow-water conditions

Higher-order FEM for nonlinear hydroelastic analysis of a floating elastic strip in shallow-water conditions VI International Conference on Computational Methods for Coupled Problems in Science and Engineering COUPLED PROBLEMS 2015 B. Schrefler, E. Oñate and M. Papadrakakis (Eds) HIGHER-ORDER FEM FOR NONLINEAR HYDROELASTIC ANALYSIS OF A FLOATING ELASTIC STRIP IN SHALLOW-WATER CONDITIONS ANGELIKI E. KARPERAKI*¹, KOSTAS A. BELIBASSAKIS*², THEODOSIOS K. PAPATHANASIOU†³ AND STILIANOS I. MARKOLEFAS‡4 *School of Naval Architecture and Marine Engineering National Technical University of Athens Heroon Polytechniou 9, Zografos 15773, Athens, Greece e-mail: ¹karperak[email protected], ²[email protected]ech.ntua.gr, http://arion.naval.ntua.gr/~kbel/ † School of Applied Mathematical and Physical Science National Technical University of Athens Heroon Polytechniou 9, Zografos 15773, Athens, Greece email: ³[email protected] ‡ Department of Mechanical Engineering Technological Educational Institute of Central Greece Psachna, Evia 34400 e-mail: [email protected] Key words: Higher-order FEM, hydroelasticity, Shallow Water Abstract. The hydroelastic response of a thin, nonlinear, elastic strip floating in shalow-water environment is studied by means of a special higher order finite element scheme. Considering non-negligible stress variation in lateral direction, the nonlinear beam model, developed by Gao, is used for the simulation of large flexural displacement. Full hydroelastic coupling between the floating strip and incident waves is assumed. The derived set of equations is intended to serve as a simplified model for tsunami impact on Very Large Floating Structures (VLFS) or ice floes. The proposed finite element method incorporates Hermite polynomials of fifth degree for the approximation of the beam deflection/upper surface elevation in the hydroelastic coupling region and 5-node Lagrange finite elements for the simulation of the velocity potential in the water region. The resulting second order ordinary differential equation system is converted into a first order one and integrated with respect to time with the Crank-Nicolson method. Two distinct cases of long wave forcing, namely an elevation pulse and an N-wave pulse, are considered. Comparisons against the respective results of the standard, linear Euler-Bernoulli floating beam model are performed and the effect of large displacement in the beam response is studied. 1 INTRODUCTION The hydroelastic interaction problem of free surface gravity waves with large, floating 1110 Angeliki E. Karperaki, Kostas A. Belibasakis, Theodosios K. Papathanasiou and Stilianos I. Markolefas. 2 bodies is found in numerous applications, ranging from the design and construction of marine structures to the response analysis of ice formations. Advances is marine technology, along with the growing need for commercial space in the condensed coastal areas has led to the rise of Very Large Floating Structures (or VLFS). Nowadays, VLFS are used, both near shore and in the open ocean, as energy plants, docking stations, storage facilities and even as floating airports and recreational amenities [1]. The analysis of such marine structures, is an important stepping stone towards robust design and construction [1-3]. Hydroelastic analysis is also relevant in the case of large ice floes under ocean wave excitation [4]. The continuous oscillatory, flexural motion undergone by large ice formations lead to their eventual splitting and disintegration. The demise of the Marginal Ice Zone (MIZ), the boundary between ice shelves and the open ocean occupied predominantly by ice floes, is linked with wave forcing [5]. In absence of a dense MIZ sums of wave energy reach the ice shelves leading to collapse events with profound environmental impact. As already mentioned, both of the above problems have their foundations set in hydroelasticity [1-4]. Their large horizontal dimensions compared to thickness, make hydroelastic effects dominant. Large floating structures are most commonly modelled in the literature as plates with zero or non-zero draft. The Kirchhoff thin plate theory is employed in the majority of works [6-8], while some consider the Reissner-Mindlin and Von Karman plate models accounting for shear deformation and large deflection effects [9-114]. For the hydrodynamic modelling, typically the linearised wave theory is utilised. When dealing with harmonic excitation eigenfunction expansion methods [12], Galerkin schemes [13] and Green functions have been employed for the solution of the hydroelastic problem in the frequency domain. However, the consideration of irregular loading dictates time domain analysis tools, such as direct integration schemes [14] and Fourier transforms [15]. Considering long wave excitation, Sturova [16], developed an eigenfunction expansion technique for the calculation of the dynamic response of a floating, thin, elastic plate of variable thickness over shallow bathymetry regions. Following the same line work, Papathanasiou et al. [17] consider a higher order finite element scheme for the solution of the transient hydroelastic problem posed by a thin, elastic, heterogeneous beam floating over shallow waters. Since large floating bodies are expected to span over great distance, the effects of variable bathymetry must also be taken into account. Belibassakis and Athanassoulis [18] derived a consistent coupled mode method for the hystroelastic analysis of a thin floating body over general bathymetry, exhibiting continuous variation. In [19] and [20] the authors extended previous work in order to account for weakly non-linear wave excitation and shear deformable bodies over general bathymetries. In the present work, the finite element method is employed for the solution of the 1D hydroelastic problem of a uniform, elastic strip floating over an uneven bottom, under shallow water conditions. The employed shallow water assumption allows for the study of tsunami impact on large floating bodies, like VLFS. The Gao beam theory [21], accounting for large defections but infinitesimal strains, is used for the approximation of the elastic strip response while the non-linear shallow water theory is chosen for the hydrodynamic model. In Section 2 of the paper the governing equations of the model in question are presented. Subsequently, in Section 3, the equivalent variational problem is derived and the proposed finite element scheme implementation is presented. In Section 4, a series of numerical results are presented 1111 Angeliki E. Karperaki, Kostas A. Belibasakis, Theodosios K. Papathanasiou and Stilianos I. Markolefas. 3 for a given configuration and variable steepness for the long wave excitation, in the form of an elevation pulse. Finally, results of the hydroelastic code are provided for the case of an Nwave excitation. 2 GOVERNING EQUATIONS In this section, the hydroelastic problem of a thin, elastic strip floating over shallow waters is presented. At the area where hydroelastic coupling is present the following equations hold,   2 2 22 4 22 ( ) (,) 2 w t rt x x x x w wt x m I D s g qxt             (1)   [() ] 0 tx x bx        , (2) where a superimposed dot denotes differentiation with respect to time. In the above equations w  , g are the water density and gravity acceleration respectively, (,)qxt denotes an external load applied to the beam and ()bx is the, possibly varying, bathymetry. Assuming that the density of the beam is e  , its Elastic Modulus E and Poisson’s coefficient  , while its thickness is  , the constants appearing in Eq. (1) are e m   the mass per width, 3 /12 re I   the rotary inertia per width and 3 11 (1 )(1 ) (1 2 ) / 12DE v v v      is the flexural rigidity per width. Finally, 21 3 (1 ) / 2sE v     is the coefficient of the nonlinear term in the Gao beam model per width [21]. Compared to the classical Euler – Bernoulli beam, the above model incorporates the effects of rotary inertia, introduced by Lord Rayleigh and the nonlinear term 22 () xx S    derived by Gao for the large deflection of a thin beam, when non-negligible stress variation in the lateral direction is considered. In addition, the pressure forcing terms   2/2 w w wx g       appearing in Eq. (1), include the nonlinear term   2 /2 wx   which is significant when the velocity x u    becomes large. In the regions where no floating structure is present, the Shallow Water Equations (SWE) are considered for the long wave propagation simulation. The equations read 0, tx x uuug      (3)   [() ] 0 tx x bx        , (4) In the following, assuming sufficient regularity and introducing the velocity potential  such that x u    , we will express equations (3) and (4) as a single evolution equation for  . For that, let us differentiate Eq. (3) with respect to t and equation (4) with respect to x , to get   2 22 2 10 2 t xt xt u ug       , (5)   22 [() ] 0 xt x bx u      , (6) 1112 Angeliki E. Karperaki, Kostas A. Belibasakis, Theodosios K. Papathanasiou and Stilianos I. Markolefas. 4 Substituting the term 2 xt   of equation (5) using (6) and setting x u    ,     2 32 2 1[() ] 0 2 xtt xt x x x g bx            (7) Integrating (7) with respect to x leads to,     2 32 2 1[() ] 0 2 xtt xt x x x g bx            (8) Select the constant () 0Ct  . Setting x u    directly in equation (3) and integrating with respect to x , we get   2 1() 2 tx x g ct        (9) Setting () 0ct  equation (9) becomes,   2 11 2 xt gg      (10) Finally, eliminating  from Eq. (8) using Eq. (10)         23 2 11 () 0 22 t t x x x x x xtx g bx            (11) Introducing the nondimensional quantities 1 x Lx   , 1/2 1/2 t gLt   , 1 L    , 1/2 3/2 gL    , where L denotes the Length of the beam, equations (1), (2) and (11) become, after dropping tildes   2 2 322 3 4 2 2 1 ( ) 2 (,) t R tx x x x t x M I K S Qxt                 , (12)   [() ] 0 tx x Bx        , (13)         23 21 1 2 () 2 0 t t x x x x x xtx Bx               , (14) where / ew M   , 12 e R w I    , (1 ) 12(1 )(1 2 ) w Ev Kv v gL    , 2 3 2(1 ) w E Sv gL   , () () bx Bx L  (,) (,) w qxt Qxt gL    , and /1L   . The bending moment and shear force inside the beam are 32 bx MK    and (15) 33 32 1 3 3( ) x Rt x x VK I S          , (16) 1113 Angeliki E. Karperaki, Kostas A. Belibasakis, Theodosios K. Papathanasiou and Stilianos I. Markolefas. 5 In equation (12), there appear four constants, namely M  , 3 R I  , 3 K  and S  . Observe that two of them scale as ()O  , while the other two are 3 ()O  . In addition, constants K and S are inversely proportional to L , the length of the beam. 2.1 Initial-Boundary Value Problem formulation In order to formulate the initial-boundary value problem of a freely-floating strip interacting with a surface wave over shallow water conditions, let us define the three nonoverlapping sets 1( ,0)   , 0 (0, )L and 2(, )L  (see Fig. 1). Figure 1: The initial boundary value problem configuration Considering the non-dimensional Eqs, (12)-(14), derived in the previous section, the IBV problem is, Find 00 :   , : ii   , 0,1, 2i , such that         23 21 1 1 1 1 1 11 2 () 2 0 t t x x x xx xt x g bx                , in 1(0, ]T (17)   2 322 3 4 2 2 0 0 0 00 2 1 00 0 () 2 (,) t R tx x x x tx MI KS Qxt                      (18) and   [() ] 0 tx x Bx        , in 0(0, ]T , (19)         23 21 1 2 2 2 2 22 2 () 2 0 t t x x x xx xt x g bx                , in 2(0, ]T . (20) with   2 1 2 i ti xi        , 1, 2i . The equations above are supplemented with the following boundary conditions, 12 (| | , ) (| | , ) 0 xx xt xt       (0, ]tT and (21) 1114 Angeliki E. Karperaki, Kostas A. Belibasakis, Theodosios K. Papathanasiou and Stilianos I. Markolefas. 6 (0,) (0,) (1,) (1,) 0 bb M tVtM tVt    , (0, ]tT . Appropriate interface conditions expressing mass and momentum conservation at the interfaces are     11 0 0 00 (0 ) (0 , ) (0 ) (0 , ) xx xx Bt Bt            , 10 00 tt xx      , (0, ]tT (22)     00 22 11 (1 ) (1 , ) (1 ) (1 , ) xx xx Bt Bt           , 02 11 tt xx      , (0, ]tT (23) The initial state for 0t describing still water conditions and zero upper surface elevation for regions 1  and 0  , while imposing an initial upper surface elevation located at a subdomain of 2  are 00 ( ,0) ( ,0) 0xx    , in 0  , (24) 11 ( ,0) ( ,0) 0 t xx   , in 1  , (25) 22 ( ,0) 0, ( ,0) ( ) t x x Gx   , in 2  , (26) 3 VARIATIONAL FORMULATION In the present section the variational formulation of the hydroelastic problem defined in Section 2.1 will be derived. Multiply Eqs. (17), (20) with 1 11 ()wH and 1 22 ()wH respectively. Multiply Eq. (18) with 2 0 ()vH (function v is not to be confused with Poisson’s ratio, a constant) and Eq. (19) with 1 00 ()wH  . Performing integration by parts it is,     00 0 22 11 11 1 1 1 1 0 00 3 1 1 1 1 1 11 1 () () 2 0 t x tx x x x x x x tx w dx w dx w B x w B x dx w dx w dx                                    , for every 1 w ,     2 3 2 3 22 0 00 00 0 1 3 3 3 32 1 3 0 000 00 2 32 1 00 0 0 000 0 0 3 () 3() 2 (,) LL L t R x tx x x L L x x x R tx x LL L L L xx t x M v dx I v dx K v dx S v dx v K I S K v v dx v dx v dx vQ x t dx                                                , for every v ,   00 0 0 0 0 0 0 0 00 [() ] () 0 LL L tx x x w dx w B x dx w B x                 , for every 0 w ,     22 22 22 2 2 2 2 3 1 2 2 2 2 22 2 () () 2 0 t x tx x L LL x x x x x tx L LL w dx w dx w B x w B x dx w dx w dx                                    , for every 2 w . 1115 Angeliki E. Karperaki, Kostas A. Belibasakis, Theodosios K. Papathanasiou and Stilianos I. Markolefas. 7 Using the boundary and interface conditions, described by Eqs. (21), (22), (23) the variational problem becomes, Find 0  and i  , 0,1, 2i , such that for every 1 () ii wH , 0,1, 2i and 2 0 ()vH it is     00 0 22 11 11 1 1 1 00 3 1 1 1 11 1 2 3 2 3 22 1 3 0 00 0 00 0 0 2 1 00 00 00 0 () 2 3 () 2 t x tx x x x x x tx LL L L t R x tx x x x x LL L tx w dx w dx w B x dx w dx w dx M v dx I v dx K v dx S v dx v dx v dx v dx w                                                          0 0 00 00 22 22 22 2 2 2 3 1 2 2 22 2 0 [() ] () 2 (,) LL tx x t x tx x x LL L L x x x tx LL dx w B x dx w dx w dx w B x dx w dx w dx vQ x t dx                                   (27) and 1 11 0 0 0 2 2 2 ( ( ,0), ) ( ( ,0), ) ( ( ,0), ) 0xw xw xw      , 0 0 0 1 11 ( ( ,0), ) ( ( ,0), ) 0 t xw xw   ,     22 2 22 ( ,0), ( ), t x w Gx w   , (,) i , 0,1, 2i being the 2 L -inner product in region i  . 3.1 Finite Element Implementation The numerical solution of the variational problem described in Eq. (27) is derived by means of the finite element method. The free water surface regions are approximated by quadratic Lagrange elements while a special element is introduced for the hydroelasticity dominated region. The reader is directed to the work of Papathanasiou et al. [17] for a more in depth analysis. The hydroelastic element incorporates fifth order Hermite polynomials for the interpolation of the beam deflection/upper surface elevation in the domain of the hydroelastic coupling and fourth order Lagrange polynomials for the interpolation of the velocity potential (see Fig. 2). The straightforward discretization of Eq. (27), and the substitution of the approximate solutions results in the following system of nonlinear ordinary differential equations, () () 0u uu uu MC K . (28) After setting uy and taking T []uyz , Eq. (28) is reduced to the first order system of nonlinear equations , () 0z zzAB . This last equation is integrated with respect to time using the Crank-Nicolson method. Figure 2: Schematic of the special hydroelastic finite element [17]. 1116 Angeliki E. Karperaki, Kostas A. Belibasakis, Theodosios K. Papathanasiou and Stilianos I. Markolefas. 8 4 NUMERICAL RESULTS Numerical results for the proposed, higher-order finite element methodology will be presented in this section. The cases of an elevated and an N-wave pulse, typical in long wave modelling, are considered. The elevated pulse is given by 2 0 00 ( ) ( )( ( ) exp( ) x w xx wxx w Gx A x x        , (29) where A is the amplitude, 0 x is the point of origin, w is the wavelength and  is a positive parameter controlling the smoothness of the initial pulse. For the isosceles N-wave profile, following Tadepalli and Synolakis [22], initial excitation is given by, 2 00 3 3 ( ) ( )sech ( ( )) 4 A Gx Adxx xx d  , (30) where d is the local depth at the origin. Figure 3, shows the corresponding initial upper surface disturbance, for each of the considered cases, that is allowed to propagate in the free surface region 2  . (a) (b) Figure 3: (a) Initial excitation in the form of an elevated pulse with 0.4Am and wavelength 265wm (b) Initial excitation in the form of an N-wave with 0.4Am and 265Lm The bathymetric profile considered in the following examples is kept flat underneath the strip, at a depth of 5 m . At a distance, equal to strip length, from the right edge of the floating body the depth is allowed to increase linearly until it is kept constant at 15 m for the rest of the 2  subdomain. Furthermore, the strip thickness is assumed uniform at 2 m , while its length is taken as 500 m . Finally, the material constants selected are material density 3 922.5 / ekg m   , water density 3 1025 / wkg m   , Young’s modulus 9 5 10E Pa  and Poisson’s ratio 0.3v . The acceleration of gravity is 2 10 / secgm . 4.1 Elevation pulse For the following analysis 100 special hydroelastic elements ( 0  ) and 10000 time steps were used for the calculation of the transient strip response. The elevation pulse parameters are 50   and 0.4A . In Figure 3, a visual representation of the upper surface elevation 1117 Angeliki E. Karperaki, Kostas A. Belibasakis, Theodosios K. Papathanasiou and Stilianos I. Markolefas. 9 solution is shown against time. At the beginning of time, the free surface disturbance set to originate in 2  , splits into two propagating waves, travelling in opposite directions. The excitation is partially reflected when it reaches the inclined seabed, and later when it impacts the strip edge (at the interface between 2  and 0  ). As the pulse reaches the strip, the hydroelastic wave begins to propagate, showing clear signs of dispersion. The waves propagating over shallower bathymetry travel at lower speeds, as expected. Figure 4: Space-time plot of the elevation pulse propagation. Mild reflection due to variable bathymetry and reflected pulses from the interaction with the floating strip are evident. Figure 5: Deflection comparison between the linear and the non-linear models, 0.2A and 265Lm . 1118