A high-order and mesh-free computational model for non-linear water waves

IS- Numerical Methods for Multiphase Flows A high-order and Mesh-free Computational Model for Non-linear Water WavesM. E. Nielsen, B. Fornberg and L. Damkilde VIII International Conference on Computational Methods in Marine Engineering MARINE 2019 R. Bensow and J. Ringsberg (Eds) A HIGH-ORDER AND MESH-FREE COMPUTATIONAL MODEL FOR NON-LINEAR WATER WAVES Morten E. Nielsen∗, Bengt Fornberg†and Lars Damkilde‡ ∗Department of Civil Engineering Aalborg University Niels Bohrs Vej 8, 6700 Esbjerg, Denmark e-mail: [email protected] †Department of Applied Mathematics University of Colorado 526 UCB, Boulder, CO 80309, USA e-mail: [email protected] ‡Department of Civil Engineering Aalborg University Thomas Manns Vej 23, 9220 Aalborg Ø, Denmark e-mail: [email protected] Key words: mesh-free methods, radial basis function-generated finite differences, high-order, node generation, water waves. Abstract. In this paper, we present the ongoing developments of a novel computational model for non-linear water waves that aims to provide a suitable framework for wave-structure interaction. The proposed model is based on radial basis function-generated finite differences, which allow for arbitrary and moving boundaries without the use of ghost nodes. In order to take advantage of the mesh-free setting, we propose a node generation strategy, suitable for moving boundaries. Numerical properties of the proposed model are investigated and finally the model is benchmarked. The proposed model is expected to provide a suitable computational framework for wave-structure interaction problems, due to its geometric flexibility and high-order nature. 1 INTRODUCTION During the latest decades, the interest in computational models for wave-structure interaction has increased, e.g. due to an increased demand for renewable energy. On the basis of this interest, several research projects [1, 2, 3] have been initiated, which seek to improve existing computational models for wave-structure interaction. Overall the improvements concern nonlinear interactions, high fidelity fluid models and computational efficiency. Ultimately, the goal is to enhance the computational models that are used for analysis and optimization of offshore structures, such as e.g. floating offshore wind turbines and wave energy converters. 1 213 Morten E. Nielsen, Bengt Fornberg and Lars Damkilde Several models have been proposed for the non-linear water wave problem, see e.g. [4, 5, 6, 7, 8, 9, 10, 11], which are based on different discretizations, ranging from traditional finite differences to spectral elements. The main concerns regarding computational modelling of non-linear water waves, in relation to wave-structure interaction, are (1) temporal stability, (2) accuracy, (3) geometric flexibility and (4) computation speed. This paper deals with the development of a novel computational model that aims to deal with the difficulties stated above, while providing a suitable framework for future coupling with structural models, e.g. the one used in [12, 13]. Specifically, we address (1) and (2) by investigating the use of radial basis function-generated finite differences (RBF-FD). The advantages of RBF-FD are geometric flexibility and high-order approximations [14, 15, 16]. Initially, the non-linear water wave problem is introduced in section 2. Next, the methodology is outlined in section 3, including a node generation strategy for the time-varying computational domain. Numerical properties of the resulting computational model are investigated in section 4, while non-linear stream functions waves are used for benchmarking of the non-linear stability and accuracy in section 5. Finally, some conclusions are stated in section 6. 2 MATHEMATICAL FORMULATION In this paper, the formulation of the fully non-linear water wave problem is based on potential flow theory and the assumption of no wave breaking to occur. In the following subsections, the governing equation and boundary conditions are presented for the two-dimensional case, which will be used in this paper. 2.1 Interior fluid domain Non-linear water waves may be modelled by means of potential flow theory, which makes it possible to describe the fluid flow by the scalar velocity potential function Φ(x,t) found by solving the Laplace equation expressed as ∇2Φ=0,x∈Ω(t),(1) where ∇=[∂ ∂x ∂ ∂z]Tis the gradient operator. Thus, equation (1) must be fulfilled within the time-varying fluid domain Ω(t), as illustrated in figure 1. 2.2 Free surface conditions For the fluid flow to evolve over time, the free surface boundary conditions must be updated over time. In the Lagrangian frame of reference, the kinematic and dynamic free surface conditions can be expressed as Dr Dt =∇Φ,x∈Γη(t),(2) DΦ Dt =−gη +1 2(∇Φ·∇Φ)=0,x∈Γη(t),(3) 2 214 Morten E. Nielsen, Bengt Fornberg and Lars Damkilde Figure 1: Fluid domain Ω(t) with specification of the free surface Γη(t) and bottom boundary Γb. where gis the gravitational acceleration, r=[xη]Tis the material node position vector, D Dt is the material derivative and Γη(t) is the subset illustrated in figure 1. 2.3 Bottom boundary condition The free-slip condition is enforced at the seabed, ∂Φ ∂n=0,x∈Γb,(4) where n=[nxnz]Tis the outward-pointing unit normal vector. 2.4 Summary of the periodic non-linear water wave problem For the sake of simplicity, we assume horizontal periodicity for all computations in this work, i.e. no lateral boundaries are present, and thus the problem is finally stated as                    ∇2Φ=0,x∈Ω(t), ∂Φ ∂z =0,x∈Γb, Dx Dt =∂Φ ∂x ,x∈Γη(t), Dη Dt =∂Φ ∂z ,x∈Γη(t), DΦ Dt =−gη +1 2(∇Φ·∇Φ) ,x∈Γη(t), (5) where also constant water depth is assumed and t= [0,T] specifies the time interval of the problem. 3 METHODOLOGY 3.1 Radial basis function-generated finite differences In radial basis function generated-finite differences (RBF-FD), the methodology is similar to the traditional finite difference method. The finite difference weights are approximated as the linear combination 3 215 Morten E. Nielsen, Bengt Fornberg and Lars Damkilde n  i=1 wiui≈Lu(x)|x=xe,(6) where Lis a linear operator, wiare the corresponding weights that approximate Lu(x) at evaluation point x=xeand nis the number of neighboring nodes used in the stencil. In this paper, the RBF used is the polyharmonic spline φ(r)=r2m−1where m∈N+. By augmenting the RBF by polynomial terms, up to an arbitrary order, several advantages may be achieved [14, 15, 16]. In order to compute the weights, we introduce an interpolant as s(x)= n  i=1 λiφ(x−xi)+ l  j=1 µjpj(x),(7) with matching conditions, that seeks to minimize far-field growth, as n  i=1 λipj(xi)=0,j=1,2, .., l, (8) where l=k+d kis the number of terms in a kth order, d-dimensional polynomial. Now, the linear system for computing the weights win Rd, when the RBF is augmented by polynomial terms up to degree k, is expressed as           p1(x1)... p l(x1) A. . ..... . . p1(xn)... p l(xn) p1(x1)... p 1(xn) . . ..... . .0 pl(x1)... p l(xn)                     w1 . . . wn wn+1 . . . wn+l           =           Lφ(x−x1)|x=xe . . . Lφ(x−xn)|x=xe Lp1(x)|x=xe . . . Lpl(x)|x=xe           ,(9) where A=Aij =φ(xi−xj) is the n×nRBF collocation matrix and 0is a l×lzero matrix. Hence, only the differentiation weights w=[w1w2... w n]Tare used for approximation of the differential operators described in section 2. In all computations conducted in this paper, the stencil sizes have been chosen to n≈2.5, and the order of the polyharmonic splines have been chosen as m=kfor kodd and m=k+ 1 for keven. 3.2 Explicit fourth-order Runge-Kutta time stepping In order to evolve the free surface variables in time, the method of lines approach is used. Specifically, the explicit fourth-order Runge-Kutta method is used, but with a slight modification. At the end of each time step, the free surface variables are reconfigured spatially by means of interpolation, such that the free surface nodes do not become too clustered along the free surface. Other strategies could have been implemented, however, only this strategy will be considered here. 4 216 Morten E. Nielsen, Bengt Fornberg and Lars Damkilde 3.3 Node generation strategy In this paper, the implemented node generation strategy is similar to the ones given in [14, 17]. However, additional steps are introduced due to the time-varying computational domain. The node generation strategy is outlined as follows: 1. Generate nodes on the boundary of the computational domain, Ω(t). 2. Generate background node set that covers Ω(t) sufficiently. 3. Remove nodes from the background node set if rs< ρ/2, where ρis a local node density measure and rsis the smallest distance to a static boundary node. 4. Perform node repelling on a subset of each of the static boundary nodes’ stencil nodes. 5. Save the repelled version of the background node set. 6. Remove nodes from the background node set if rd< ρ/2, where ρis a local node density measure and rdis the smallest distance to a dynamic boundary node. 7. Perform node repelling on a subset of each of the dynamic boundary nodes’ stencil nodes. 8. Remove nodes from the background node set if they are located outside Ω(t). 9. Repeat step 6-8 at each substep, where the saved background node set from step 5 is re-used each time. In the strategy outlined above, static and dynamic boundary nodes refer to their variation over time. Hence, only nodes near moving boundaries are updated at each substep. In this paper, the node generation algorithm from [17] is used to generate the background node set. 3.4 Solution procedure The overall solution procedure is outlined in algorithm 1. Using an explicit time stepping method, e.g. the RK4 as presented previously, the time derivatives need to be computed, as outlined in section 2 and 3. First the initial- and boundary conditions Φ0 (x,η)0,η0 x0,x 0are provided as input. These are used to solve the Laplace equation in (1) for the velocity potential Φ(x, z, t)=Φ t (x,z)within Ω(t). Next, the particle velocities are computed at the free surface, which are then used to compute the temporal derivatives given in equations (2,3). In this manner, the free surface variables are stepped forward in time. Before proceeding, the new boundary conditions, Φt (xt,ηt xt),ηt xt,x t, are interpolated back to their initial horizontal position, x0, and the next time step is initialized from this spatial configuration, i.e. Φt (x0,ηt x0),ηt x0,x 0. 5 217 Morten E. Nielsen, Bengt Fornberg and Lars Damkilde Algorithm 1: Overall solution procedure Input: Φ0 (x,ηx0)0,η0 x0,x 0 1t:= 0 2while t<T do 3Φt+∆t (x,η)t+∆t,ηt+∆t xt+∆t,x t+∆t←RK4(L,Φt (x0,ηt x0),ηt x0,x 0); % Time stepping 4ηt+∆t x0←interp(ηt+∆t xt+∆t,x t+∆t,x 0); % Interpolate free surface elevation 5Φt+∆t (x0,ηt+∆t x0)←interp(Φt+∆t (x,η)t+∆t,x t+∆t,x 0); % Interpolate velocity potential 6t:= t+∆t 7end Output: Φ0 (x,η)0,···,ΦT (x0,ηT),η0 x0,···,ηT x0,x 0 4 NUMERICAL ACCURACY AND STABILITY First, the accuracy of computing the vertical particle velocity is examined for different ranges of non-linearities. Next, the temporal stability is investigated and it is demonstrated how to improve the stability by means of a previously proposed method [18]. 4.1 Accuracy The spatial accuracy is illustrated by the relative error between vertical particle velocities, at the free surface, computed by stream function theory and the proposed model, respectively. The wave steepness (H/L), node density (ρ) and order of augmented polynomials (P) are varied, while the wave length L=2πis kept constant. The results are illustrated in figure 2 along with an example of the node set for the corresponding wave steepness. From figure 2, it is noticed that the convergence rate decreases as the wave steepness approaches the breaking limit. 4.2 Linear stability analysis The purpose of the linear stability analysis is to investigate the temporal stability when different polynomials are used, and to test if instabilities can be removed. First, we express the linearized and semi-discrete counterpart of the system in (5) as D Dt   x η Φ =J  x η Φ =  00Lx 00Lz 0−gI0   x η Φ ,(10) where Iis the identity matrix and Lzis derived from the relation DzΦ|x=r= [(Dz)11Φ1+(Dz)12Φ2]=(Dz)11 −(Dz)12L−1 22 L21Φ1=LzΦ1,(11) where Lis the differential operator for the Laplace problem, introduced in section 2, Dzis the differentiation matrix with respect to the vertical direction and the matrices are decomposed as L=L11 L12 L21 L22 ,(12) 6 218 Morten E. Nielsen, Bengt Fornberg and Lars Damkilde kh =2π,H/L =0.014 0.1 0.15 0.2 0.25 0.3 0.35 0.4 10-6 10-4 10-2 100 kh =2π,H/L =0.07 0.1 0.15 0.2 0.25 0.3 0.35 0.4 10-6 10-4 10-2 100 kh =2π,H/L =0.126 0.1 0.15 0.2 0.25 0.3 0.35 0.4 10-6 10-4 10-2 100 Figure 2: Convergence plots of the relative error computing ∂Φ ∂z at the free surface for different H/L ratios and augmented polynomials, while using a node density approximately equal to ρ. where the subscript indices define whether the nodes are on the free surface (1) or not (2). Lx is computed in a similar manner as in equation (11). Eigenspectra of the system in equation (10) are illustrated in figure 3.a, where a quasi-uniform node set has been used and also different orders of augmented polynomials. As noticed, the eigenvalues in figure 3.a has positive real parts, which leads to temporal instability. However, as described in e.g. [18], convective PDEs 7 219 Morten E. Nielsen, Bengt Fornberg and Lars Damkilde may be stabilized by adding a slight amount of artificial dissipation, known as hyperviscosity. In this work, a third order Laplacian is used along with a tuning parameter γ, which depends on the discretization used. The effect of adding hyperviscosity to the free surface elevation and velocity potential, in equation (10), is illustrated in figure 3.b. Here γ=6·10−6has been used for all cases, although a more suitable choice exist for each individual discretization. a) without hyperviscosity -0.2 -0.1 0 0.1 0.2 -10 -5 0 5 10 b) with hyperviscosity -10 -8 -6 -4 -2 0 2 -10 -5 0 5 10 -2 -1 0 1 2 -10 -5 0 5 10 -10 -8 -6 -4 -2 0 2 -10 -5 0 5 10 -4 -2 0 2 4 -10 -5 0 5 10 -10 -8 -6 -4 -2 0 2 -10 -5 0 5 10 -0.2 -0.1 0 0.1 0.2 -10 -5 0 5 10 -10 -8 -6 -4 -2 0 2 -10 -5 0 5 10 Figure 3: Eigenspectra of Jin equation (10) using a quasi-uniform node set and augmented by up to the Pth order polynomial. 8 220 Morten E. Nielsen, Bengt Fornberg and Lars Damkilde 5 RESULTS Stream function waves of varying steepness are used as benchmark for the non-linear stability and accuracy. Each case is simulated for five periods with the time step ∆t=0.025s≈T/80, where Tis the wave period. The maximum relative errors between the stream function solution and proposed model are shown in figure 4. For each discretization, a suitable γ-value is determined by linear stability analyses. kh =2π,H/L =0.014 0.1 0.2 0.3 0.4 10-4 10-3 10-2 10-1 100 kh =2π,H/L =0.07 0.1 0.2 0.3 0.4 10-4 10-3 10-2 10-1 100 kh =2π,H/L =0.126 0.1 0.2 0.3 0.4 10-4 10-3 10-2 10-1 100 Figure 4: Convergence plots showing the relative error of the free surface elevation after five periods when using a 7th order polyharmonic spline augmented by up to 6th order polynomials. Once again, it is noticed that the convergence rate tend to decrease as the wave steepness increase. As expected, it drops to the same convergence rate that was noticed in figure 2. However, as quasi-uniform nodes were used, it can be expected that node refinement near the free surface will improve the accuracy without increasing the total number of nodes used. 6 CONCLUSIONS A novel computational model for non-linear water waves is introduced and the methodology outlined. The main novelty is the application of RBF-FD, which provide high-order approximations along with geometric flexibility for the moving free surface. Additionally, a modification of a previously developed node generation strategy is introduced, such that the computational domain is updated as the boundaries evolve over time. Finally, it is shown that stability can be achieved by adding hyperviscosity to the problem without influencing the convergence rate. REFERENCES [1] DeRisk, http://www.derisk.dk, [Online; accessed 11-February-2019] (2015). [2] MIDWEST, https://project.inria.fr/midwest/, [Online; accessed 11-February-2019] (2016). [3] CCP-WSI, https://www.ccp-wsi.ac.uk/, [Online; accessed 11-February-2019] (2016). 9 221