- Email: [email protected]
LRBF-DQ
Ocean Modelling 65 (2013) 1–10
Contents lists available at SciVerse ScienceDirect
Ocean Modelling journal homepage: www.elsevier.com/locate/ocemod
Numerical modeling of long waves in shallow water using LRBF-DQ and hybrid DQ/LRBF-DQ A. Khoshfetrat 1, M.J. Abedini ⇑ Dept. of Civil and Environmental Engineering, Shiraz University, Shiraz, Iran
a r t i c l e
i n f o
Article history: Received 13 August 2012 Received in revised form 12 January 2013 Accepted 20 January 2013 Available online 15 February 2013 Keywords: Differential quadrature Local RBF-based differential quadrature Shallow water Oresund strait Hybrid methods
a b s t r a c t A computational algorithm based on coupling the differential quadrature (DQ) method and the local radial basis function-based differential quadrature (LRBF-DQ) method is used to solve the shallow water equations. The proposed method that is named ‘hybrid DQ/LRBF-DQ’ has the advantages of DQ such as simplicity and low computational cost and the advantages of LRBF-DQ, like mesh free and Kronecker delta function properties. Contrary to the DQ method, this method can be used in problems with irregular domain with less computational cost and less sensitivity to shape parameter(s) compared to the LRBF-DQ method. The proposed method is employed to simulate the complex and irregularly shaped waterway of the Oresund strait which is located between Sweden and Denmark. Implementation of the proposed method leads to good results for water depth. In comparison with the LRBF-DQ method, the results show that the proposed method is less sensitive to the shape parameter and the range of acceptable shape parameter in the proposed method is wider. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction The calculation of tidal elevations and currents in coastal regions constitute an important part of any coastal management study. As solving the governing partial differential equations (PDEs) analytically is impossible, an approximate numerical solution must thus be found using a suitable numerical method. The most popular methods are the finite difference (Rasulov et al., 2005; Skiba and Filatov, 2008; Wang and Liu, 2011), the finite element (Hanert et al., 2005; Cotter et al., 2009; Kesserwani and Liang, 2010), the finite volume (Castro et al., 2008; Daoud et al., 2008; Adamy et al., 2010) and the method of characteristics (Toda et al., 2009). The Differential Quadrature Method (DQM) proposed by Bellman and Casti (1971), Bellman et al. (1972) is an alternative numerical method for the numerical solutions of linear and nonlinear partial differential equations. Using concepts borrowed from polynomial approximation and linear vector space, Shu and Richards (1992) presented the polynomial-based differential quadrature (PDQ), in which the weighting coefficients of any order derivative are calculated by a simple algebraic formulation or a recurrence relationship without any restriction on the choice of grid distribution. Later, Shu and Chew (1997) also developed the Fourier-expansion-based differential quadrature (FDQ) where the ⇑ Corresponding author. Tel.: +98 711 6474604; fax: +98 711 6473161. 1
E-mail address: [email protected] (M.J. Abedini). Current address: Khorasgan Branch, Islamic Azad University, Isfahan, Iran.
1463-5003/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ocemod.2013.01.006
function is approximated by the Fourier series expansion. The details of PDQ and FDQ as well as their applications in different disciplines can be found in Shu (2000). Due to its efficiency and accuracy, DQM has been widely employed in many areas of industry and mathematics (Wu and Shen, 2004; Wu and Kong, 2005; Wu and Liu, 2005; Ren and Wu, 2006). Gutierrez et al. (1994) proposed the implementation of DQM in the field of ocean engineering. Hashemi et al. (2006, 2007, and 2008) used DQM as a rapid and accurate method for numerical simulation of one dimensional shallow water long waves. They showed the advantages for DQM compared with other conventional methods, namely, the accuracy, simplicity and minimal computational effort. For the implementation to become a practical alternative to existing integration methods, further research is necessary. As an example, the method should address the extension of the method to a tidal problem in an irregular domain. On the other hand, the use of mesh free (MFree) methods has been a subject of research interest in physical and engineering problems. In the MFree methods, the field function is constructed entirely from a set of arbitrarily distributed nodes (Belytschko et al., 1996; Gu, 2005). The nodal connectivity is imposed with the use of influence domains instead of elements being conventional in the finite element method (Zienkiewicz, 1989). These methods have been recently used to model shallow water equations (Du, 2000; Rodriguez-Paz and Bonet, 2005; Darbani et al., 2010; Hon et al., 1999; Wong et al., 2002). Homayoon et al. (2013) used RBF-DQ that is an MFree method, for modeling the shallow water equations. This method suffers
2
A. Khoshfetrat, M.J. Abedini / Ocean Modelling 65 (2013) 1–10
from two main drawbacks which are high computational cost and sensitivity to shape parameter. Shu et al. (2003) introduced Local RBF-DQ (LRBF-DQ) in which the two drawbacks of RBF-DQ are reduced. In this paper, LRBF-DQ method will be used to model shallow water equations in Oresund strait. It is often desirable and beneficial to combine two established numerical methods to exploit their advantages while avoiding their disadvantages (Gu and Liu, 2005). Khoshfetrat and Abedini (2012) introduced a new method, which is named hybrid DQ/LRBF-DQ, through coupling the conventional DQ and the LRBF-DQ to solve PDEs. For this purpose, the computational domain is divided into a few rectangular shapes and some irregular shapes. In such a domain decomposition process, a high percentage of the computational domain will be covered by regular shapes thus taking advantage of conventional mesh-based DQM eliminating the need to implement local RBF-DQ over the entire domain but only on a portion of it. By this method, we have the advantages of DQ like simplicity, high accuracy, and low computational cost and the advantages of LMQRBF-DQ like mesh free and Kronecker delta function properties. While irregular geometry has been the main impetus to use domain decomposition in this study, the multi-domain technique has been implemented in other studies for totally different purposes (Shu et al., 1998; Zong et al., 2005). Shu et al. (1998) used multi-domain technique to remove the effect of singularity and Zong et al. (2005) proposed domain decomposition to address material discontinuity in the computational domain. It is expected that the computational cost in the proposed methodology be less than the LRBF-DQ method. It will also be shown that the sensitivity to shape parameter in the proposed methodology is less than the LRBF-DQ method. We will use the hybrid DQ/LRBF-DQ method to model the shallow water equations in Oresund strait and the results will be compared with observed data, the output of MIKE21 model and that obtained by LRBF-DQ method. 2. Shallow water equations (SWEs) Shallow water free surface flows are usually described by the SWEs which can be derived from the depth averaged Navier– Stokes equations over the flow depth. Assuming the fluid density is constant and the pressure distribution is hydrostatic, the continuity equation can be written as:
@h @P @Q þ þ ¼0 @t @x @y
ð1Þ
u¼
1 h
Z
zs
Udz;
v¼
zf
1 h
Z
zs
Vdz
ð5Þ
zf
where U and V are the local velocity components and Zf and Zs are the bed and surface elevations. In Eq. (3) g is the acceleration of gravity; n is the Manning roughness coefficient; f = 2xSink is the Coriolis parameter; x is the angular velocity of the earth rotation; k is the latitude of the studied area; q is the water density; qair is the air density; Cw is the wind friction factor; Uw and Vw are the wind speed components in the x and y directions. The governing equations (1)–(3) are subject to appropriate initial and boundary conditions. Usually, the domain is limited by the bottom, the water surface, the land boundaries and the boundaries which could be purely imaginary that is named open boundaries. The boundary conditions on bottom and free surface have been considered in deriving the SWEs. Strictly speaking, the boundary condition for velocity on a land boundary is a no-slip condition: ~ ¼ 0, in which U ~ ¼ ðU; VÞ is the velocity vector. However, due to U the turbulence and the existence of a boundary layer, the velocity near a wall quickly becomes non-zero, and the no-slip condition is often replaced by the impermeability or free slip boundary condi~~ tion: U; n ¼ 0, where ~ n is the unit outward normal vector on the land boundary. On the open boundaries, additional information is required on the sea surface elevation as a function of time: Z s ¼ Z s ðtÞ, where Z⁄ is the specified sea surface elevation on the ~ ¼ ð0; 0Þ, which correopen boundary. The initial conditions are U sponds to the still water condition and Z s ¼ cte, that is the mean sea level in the computational domain. 3. Hybrid DQ/LRBF-DQ method 3.1. Description In the hybrid DQ/LRBF-DQ method, the domain of problem is decomposed into some regular and irregular subdomains. Consider a typical two-dimensional problem as shown in Fig. 1a. The problem domain consists of five subdomains X1 ; X2 ; X3 ; X4 ; X5 . If the point (xi, yj) is located inside region X1, then the conventional DQM will be used to approximate the derivative with respect to x, ux and the derivative with respect to y, uy at (xi, yj) via:
ux ðxi ; yj Þ ¼
NX X wxi;k uðxk ; yj Þ
ð6aÞ
k¼1
uy ðxi ; yj Þ ¼
NY X wyj;k uðxi ; yk Þ
ð6bÞ
k¼1
The momentum equation in x-direction can be expressed as: 2
!
@Z f @P P @h 2P @P PQ @h Q @P P @Q þ gh 2 þ þ þ þ gh @t h @x h2 @y h @y h @y @x h @x qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi gn2 P P2 þ Q 2 q þ fQ air C w U w U 2w þ V 2w ¼ 0 7=3 q h
ð2Þ
where NX, NY are, respectively, the number of mesh points in region X1 in the x and y directions, wxi,k, wyj,k are the DQ weighting coefficients in the x and y directions. In order to determine the above weighting coefficients, many kinds of test functions have been used including polynomials, the sine functions, the cosine functions and various orthogonal polynomials. The Lagrange interpolation function is widely used for this purpose, namely
For the y-direction, the momentum equation is:
lk ðxÞ ¼
2
@Z f @Q PQ @h Q @P P @Q Q @h 2Q @Q þ 2 þ þ þ ðgh 2 Þ þ þ gh @t @y h @y h @y @y h @y @y h h qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi gn2 Q P2 þ Q 2 q þ fP air C w V w U 2w þ V 2w ¼ 0 þ 7=3 q h
ð3Þ
where t is time; h is the flow depth; P and Q are the discharge per unit width in x and y directions that can be expressed as:
P ¼ uh;
Q ¼ vh
ð7Þ
Substituting Eq. (7) into Eq. (6a), we can get the following formulas to directly calculate the weighting coefficients of the firstand second-order derivatives, i.e.
wxi;j ¼
ð4Þ
In which, u and v are the depth averaged velocities in x and y directions which are defined respectively as:
NX Y ðx xi Þ x xi i¼1;i–k k
wxi;i ¼
NX Y 1 xi xk ; for j–i xj xi k¼1;k–i;j xj xk NX X
1 x xk k¼1;k–i i
ð8Þ
ð9Þ
A. Khoshfetrat, M.J. Abedini / Ocean Modelling 65 (2013) 1–10
(a)
3
(b)
Fig. 1. Domain decomposition for Hybrid DQ/RBF-DQ method.
The weighting coefficients in the y direction can be obtained in a similar way. If the point ðxl ; yl Þ is located in the irregular subdomains X2 ; X3 ; X4 ; X5 (containing the interface regions), local RBF-DQ is used to find the derivative with respect to x, ux and the derivative with respect to y, uy at ðxl ; yl Þ by:
ux ðxl ; yl Þ ¼
N X wxl;m uðxm ; ym Þ
ð10aÞ
m¼1
uy ðxl ; yl Þ ¼
N X wyl;m uðxm ; ym Þ
ð10bÞ
m¼1
where N is the number of supporting points in the neighborhood of ðxl ; yl Þ; wxl;m and wyl;m are the LRBF-DQ weighting coefficients and u(xm,ym) is the function value at all supporting points (xm,ym). The weighting coefficients in Eq. (10) are determined by using the multi-quadric (MQ) RBFs as test functions. MQ RBFs can be expressed as follows:
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi up ðXÞ ¼ uðjjX X p jjÞ ¼ uðrÞ ¼ r2 þ c2
ð11Þ
where X ¼ ðx; yÞ; r ¼ jjX X p jj is the Euclidean norm, and c is a free shape parameter. Substituting the set of MQ radial basis functions into Eqs. (10), one can obtain: N @ up ðxl ; yl Þ X ¼ wxl;m up ðxm ; ym Þ @x m¼1
ð12aÞ
N @ up ðxl ; yl Þ X wyl;m up ðxm ; ym Þ ¼ @y m¼1
ð12bÞ
For p ¼ 1; 2; . . . ; N. Equation systems (12a) and (12b) are used to compute the weighting coefficients wxl;m and wyl;m , respectively. For the case when the total number of points in the computational domain are used in Eqs. (10) and (12), then the derivatives is said to be approximated by the RBF-DQ method. Now, we can combine the DQ and LRBF-DQ equations used in hybrid DQ/LRBF-DQ method, into one system of equations:
ux ðxi ; yi Þ ¼
Nix X wxi;j uðxj ; yj Þ
ð13aÞ
j¼1
uy ðxi ; yi Þ ¼
Niy X wyi;j uðxk ; yk Þ
ð13bÞ
k¼1
where (xi,yi) is the point at which the derivatives are to be approximated. These are the formulas for derivative approximation in the proposed hybrid DQ/LRBF-DQ method. To determine the Nix, Niy and points (xj,yj), (xk,yk) in equations (13a, b), If the point (xi,yi) is located inside region X1, then the conventional DQM will be used to
compute the weighting coefficients in Eqs. (13a, b), so the Nix and Niy are, respectively, the number of rectangular mesh points (including the boundary points) in the x and y directions (Fig. 1b), and yj = yi, xk = xi in Eqs. (13a) and (13b), respectively. If the point (xi,yi) is located in the irregular subdomains X2,X3,X4,X5 (containing the interface regions), the Local RBF-DQ is used to find the weighting coefficients in Eqs. (13a) and (13b). Nix = Niy are the number of supporting points surrounding the point under consideration and uðxj ; yj Þ ¼ uðxk ; yk Þ is the function value at the supporting point. Supporting points may be located either in regular or irregular regions. Once again, spatial and temporal connectivity of state variable(s) between regular and irregular domains will be taken care of by points common to both regions. By using this method, we have the low cost DQM in a large part of domain and a mesh free Local RBF-DQ in the small regions near the boundaries for implementation of non-rectangular geometry which cannot be considered by the conventional DQM (Fig. 1b). 3.2. Advantages over traditional methods - The proposed hybrid approach and its variants (i.e., DQ and RBF-DQ) are categorized as higher order accurate due to reservation of a large number of terms in the Taylor series expansion. The exact order depends on the number of nodes considered to expand the derivative. - For a conventional triangular linear element extensively used in FEM and FVM, the state variable(s) would vary linearly between nodal points, while in RBF-DQ with multi-quadric basis function, the state variable(s) would vary in an exponential manner between nodal points. The major implication of variation in state variable(s) between nodal points concerns with rate of convergence. Multi-quadric RBFs have an exponential rate of convergence, while the conventional FEM and FVM implementation has linear rate of convergence. That is why it is argued that with much fewer nodes; the results obtained compared well with analytical and benchmark solutions of other numerical schemes. In other words, with a similar node distribution in RBF-DQ and FE, by similarly refining the nodes, the RBF-DQ results converge to the correct solution more rapidly than those of FEM or FVM. Furthermore, it might be interesting to know that one more advantage of RBF-based method as summarized by Wu (1992) is that it is space dimension independent in the sense that the convergence order is of O(hd + 1) where h is the density of the collocation points and d is the spatial dimension. - It is a well-known fact that traditional numerical techniques such as FDM, FEM, and FVM are strongly dependent on mesh properties. In light of this, mesh generation takes up very much time for complex geometries while resorting to these methods for obtaining numerical solution of continuum mechanics problems. Indeed, the mesh generation is considered to be the most expensive part of the simulation in these methods. Mesh-free methods (also called mesh-less methods) such as RBF-DQ have
4
A. Khoshfetrat, M.J. Abedini / Ocean Modelling 65 (2013) 1–10
generated considerable interest in recent years due to the need to overcome the high cost of mesh generation associated with the human labor in mesh-based methods. In mesh-free methods, such mesh generation is not required and instead, node generation is done which is much easier. In node generation, node connectivity and the area between the nodes is not important as opposed to in mesh generation. As an example, in FEM, if a triangular element becomes so narrow, poor results is obtained. However, in RBF-DQ, the characteristics of the area between the nodes are of no importance. As a result, node generation in RBF-DQ is very easier than mesh generation in FEM. Actually, the notion ‘‘mesh-free’’ is more suitable than ‘‘meshless’’, as in these methods, the nodes could have more degree of freedom compared to mesh-based methods. In short, the cost of mesh generation is a big issue in mesh-based methods. Generally speaking, in contrast to mesh-based methods where the grid resolution and/or choice of basis function are the primary determinant of accuracy, the shape parameter plays an important role in the accuracy of RBF-based methods in addition to the node distribution/density. The salient feature of our proposed hybrid numerical scheme is to diminish the impact of shape parameter on numerical accuracy. - In mesh-based methods such as FEM, the interpolating function will be defined over a typical element, while in mesh-free methods such as RBF-DQ; the interpolating function will be defined over a typical node as opposed to a typical element. The major implication of this change in focus from ‘‘element to node’’ relates to node refinement. In mesh-free methods, node refinement can be handled very easily, while in meshbased methods, re-meshing would be quite cumbersome for irregular domains if the numerical exercise requires doing so. - In FEM (or FVM), definition of interpolating function over a typical element is highly localized, while in RBF-DQ, its definition has a global nature. Indeed, the state variable has nonzero values over a typical element in FEM (or FVM) and zero value elsewhere. This won’t be the case in RBF- DQ while using infinitely smooth RBFs such as multi-quadric. The shape parameters are intended to control the impact of neighboring nodal points.
4. Numerical formulation of SWEs In this paper, the hybrid DQ/LRBF-DQ and LRBF-DQ methods are adapted to solve the shallow water equations (1)–(3) subject to the initial and boundary conditions as described in Section 2. Eqs. (1)– (3) can be approximated as: tþ1
hi
t
Nix X hi X wxij Ptj þ wyik Q tk ¼ 0 þ Dt j¼1 k¼1 Niy
ð14Þ
0 !2 1 Niy Nix Nix P tþ1 P ti @ t P ti AX 2P t X Pt Q t X t t i wxij hj þ t i wxij P tj i t 2i wyik hk þ ghi t Dt hi hi j¼1 ðhi Þ k¼1 j¼1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 t Niy Niy tX Nix X gn P ðP ti Þ þ ðQ ti Þ Q ti X P i t þ t wyik P tK þ ti wyik Q tk þ ghi wxij Z fj þ t 7=3 hi k¼1 hi k¼1 ðhi Þ j¼1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q t 2 2 ð15Þ fQ i þ air C w U tw ðU tw Þ þ ðV tw Þ ¼ 0
q
Nix Nix Nix Q itþ1 Q ti Pti Q ti X Qt X Pt X t wxij hj þ ti wxij Ptj þ it wxij Q tj t 2 Dt hi j¼1 hi j¼1 ðhi Þ j¼1 0 !2 1 Niy Niy Niy X Q ti AX 2Q t X t t t @ þ ghi wyik hk þ t i wyik Q tk þ ghi wyik Z fk t hi hi k¼1 k¼1 k¼1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi gn2 Q ti ðPti Þ2 þ ðQ ti Þ2 q t þ þ fPi air C w V tw ðU tw Þ2 þ ðV tw Þ2 ¼ 0 t 7=3 q ðhi Þ
ð16Þ where wxij and wyik are the hybrid DQ/LRBF-DQ or the LRBF-DQ weighting coefficients that are determined according to previous section. Nix and Niy are the number of supporting points in the LRBF-DQ method or the number of points in x or y direction in DQ method. The subscript and superscript of h, P and Q denote the location and time, in which the equations are approximated, respectively. It can be seen that the equations are discretized by explicit method in time, or the state variables in time step t + 1 are unknown.
Fig. 2. Location of Oresund strait (inside the circle).
A. Khoshfetrat, M.J. Abedini / Ocean Modelling 65 (2013) 1–10
~~ For satisfying the impermeability condition: U: n ¼ 0 on the land boundary, the flow velocity on this boundary needs to be adjusted so that the normal flow velocity is equal to zero. Here, we have the discharge per unit width instead of flow velocity. In each time step, the computed discharges on the land boundaries (P and Q) are updated as:
P ¼ Pn2y Qnx ny
ð17Þ
Q ¼ Qn2x Pnx ny
ð18Þ
where nx and ny are the x and y components of the normal vector on the land boundary, respectively. Eqs. (17) and (18) eliminate the normal component of the calculated discharge on the land boundary.
5
5. Solving SWEs in Oresund strait Oresund strait, located between Denmark and Sweden, connects Baltic Sea to Kattegat gulf (Fig. 2). In 1994, the construction of a fixed link between Copenhagen (Denmark) and Malmo (Sweden) as a combined tunnel, bridge and reclamation project commenced. Severe environmental constraints were enforced to ensure that the environment of the Baltic Sea remains unaffected by the link. To meet the environmental constraints and to monitor the construction work, a major monitoring program was set up. The resulting data were used to calibrate the DHI’s three-dimensional model, MIKE 3 and its two-dimensional model, MIKE 21. MIKE21 Hydrodynamic Software considers Alternating Direction Implicit (ADI) technique to integrate the equations for mass and momentum conservation in the space–time domain (DHI, 2005).
Fig. 3. The limits of modeling domain (inside the dashed rectangular).
6
A. Khoshfetrat, M.J. Abedini / Ocean Modelling 65 (2013) 1–10
Here, the data reported in MIKE 21 are used to evaluate the LRBFDQ and hybrid DQ/LRBF-DQ methods and the modeling results are compared with MIKE 21 results. The part of Oresund strait located in the 63900m 83700m dashed lines rectangular that is shown in Fig. 3 is modeled. The x and y axis are selected in a way that at open boundaries, the flow
direction is perpendicular to boundary of domain. The origin is located at Easting 337100 and Northing 6122900. The domain is discretized as shown in Fig. 4. The grid spacing in x and y directions is 900m. These grid points which are considered in MIKE 21 model, is used for modeling the Oresund strait by the LRBF-DQ and the hybrid DQ/LRBF-DQ methods.
Fig. 4. Grid points considered in modeling domain and the location of measuring stations.
Fig. 5. Comparison of LRBF-DQ and measured and MIKE21 simulated results at Drogden station.
A. Khoshfetrat, M.J. Abedini / Ocean Modelling 65 (2013) 1–10
In addition to the bathymetry of all grid points that is available, the time series of the following data has been measured in 1800 s time intervals in a 12 days period from 12/2/1993 to 12/ 13/1993:
7
- Water depth at two points of each open boundary (WL1, WL2, WL3 and WL4) in Fig. 4. These data are used to compute the water depth at all grid points in the open boundaries by linear interpolating method.
Fig. 6. Comparison of LRBF-DQ and measured and MIKE21 simulated results at NdrRoese station.
Fig. 7. Domain decomposition of Oresund strait for hybrid DQ/LRBF-DQ method.
8
A. Khoshfetrat, M.J. Abedini / Ocean Modelling 65 (2013) 1–10
- Water depth at Drogden and NdrRoese stations (Fig. 4) that are used for model validation. - Wind speed and direction at strait. The water elevation and discharge per unit width at t = 0 are considered 0.38m and zero, respectively (initial conditions). The LRBF-DQ method has been used to model SWEs in Oresund strait and the results are shown in Figs 5 and 6. In these figures, the computed water elevations at Drogden and NdrRoese stations during the 12 days measuring period are plotted. To evaluate the LRBF-DQ method in modeling SWEs, the measured water elevations and the results of MIKE 21 model are plotted in these figures as well. The time step of modeling is 100 s and the number of supporting points used to compute the shape functions is 200. The various values have been considered for the shape parameter and the best results are shown in Figs. 5 and 6. Implementation of the hybrid DQ/LRBF-DQ method involves the decomposition of computational domain into some regular and irregular subdomain. This decomposition for Oresund strait is shown in Fig. 7. It can be seen that each regular region contains four nodes in each of the x and y directions. It means that in each regular subdomain, the derivatives in four interior nodes are approximated by the low-priced DQ method. According to domain
decomposition shown in Fig. 10, in 34% of the nodes, the DQ method will be used for numerical modeling when using hybrid DQ/ LRBF-DQ method. This is why the experience on SWEs modeling in Oresund strait show that if the hybrid DQ/LRBF-DQ method is used, the computational time will be reduced to 84% of the computational time involved when the LRBF-DQ method is used. The computed time series of water depth by the hybrid DQ/LRBF-DQ method at Drogden and NdrRoese stations are plotted in Figs 8 and 9, respectively. The performance of hybrid method can be evaluated by comparing the numerical results, obtained from this method, with actual measurements and MIKE21 results. Here, in hybrid DQ/LRBF-DQ method, the modeling time interval is considered as 100 s and the number of supporting points is 200. The accuracy of results depends on shape parameter value selected. The various values have been considered for the shape parameter and the best results are plotted in Figs. 8 and 9. To assess and evaluate the performance accuracy of various approaches, the following error norm is considered:
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N N X X 1u Þ2 = u 2 Error ¼ t ðu u N i¼1 i¼1
Fig. 8. Comparison of hybrid DQ/LRBF-DQ and measured and MIKE21 simulated results at Drogden station.
Fig. 9. Comparison of hybrid DQ/LRBF-DQ and measured and MIKE21 simulated results at NdrRoese station.
ð19Þ
A. Khoshfetrat, M.J. Abedini / Ocean Modelling 65 (2013) 1–10
9
model and good agreement is achieved. The proposed method has the advantages cited in Section 3.2 compared to the MIKE21 model in which a traditional numerical method is used. The results show that the sensitivity of results to the shape parameter is reduced when using hybrid DQ/LRBF-DQ method. Additionally, the computational time is reduced in hybrid DQ/LRBF-DQ method. Results of this research show that the hybrid DQ/LRBF-DQ method is an attractive alternative for DQ and LRBF-DQ methods in ocean modeling problems that reduces their shortcomings.
References
Fig. 10. Variation of error with shape parameter (c) for various methods at Drogden station.
Fig. 11. Variation of error with shape parameter (c) for various methods at Drogden station.
where N is the number of data measured at stations, u is the simulated value of state variable obtained by a particular numerical is the measured data . method, and u To evaluate the influence of the value of shape parameter on accuracy of results obtained by the LRBF-DQ and the hybrid DQ/ LRBF-DQ methods, the variation of error with various values of shape parameter at Drogden and NdrRoese stations are plotted in Figs. 10 and 11, respectively. These figures show that the hybrid DQ/LRBF-DQ method is less sensitive to the shape parameter, for the reason that range of acceptable shape parameters in which the solution is convergent and steady, is larger and the error vs. c curve is less fluctuating.
6. Conclusions In this paper, the hybrid DQ/LRBF-DQ method is used in ocean modeling for the first time. This method prepares the opportunity of using the simple, accurate and cheap DQ method for problems with irregular domains, while the conventional DQ method can be used only in regular domains. Although the LRBF-DQ method is able to model problems with such irregular domains but it involves performing a large amount of computations. Another inconvenience in using LRBF-DQ is its sensitivity to the value selected for shape parameter. To illustrate the advantages of the proposed method, the SWEs are simulated in Oresund strait. The computed results by the LRBFDQ and the hybrid DQ/LRBF-DQ methods are compared with actual measurements and good agreement is indicated. The computed results are also compared with the results obtained by the MIKE21
Adamy, K., Bousquet, A., Faure, S., Laminie, J., Temam, R., 2010. A multilevel method for finite volume discretization of the two-dimensional nonlinear shallowwater equations. Ocean Model. 33(3–4), 235–256. Bellman, R., Casti, J., 1971. Differential quadrature and long term integration. J. Math. Anal. Appl. 34, 235–238. Bellman, R., Kashef, B.G., Casti, J., 1972. Differential quadrature: a technique for the rapid solution of nonlinear partial differential equations. J. Comput. Phys. 10, 40–52. Belytschko, T., Krongauz, Y., Organ, D., Fleming, M., Krysl, P., 1996. Meshless methods: an overview and recent developments. Comput. Methods Appl. Mech. Eng. 139 (1), 3–47. Castro, M.J., García-Rodríguez, J.A., González-Vida, J.M., Parés, C., 2008. Solving shallow-water systems in 2D domains using finite volume methods and multimedia SSE instructions. J. Comput. Appl. Math. 221 (1), 16–32. Cotter, C.J., Ham, D.A., Pain, C.C., 2009. A mixed discontinuous/continuous finite element pair for shallow-water ocean modelling. Ocean Model. 26 (1–2), 86–90. Darbani, M., Ouahsine, A., Villon, P., Naceur, H., Smaoui, H., 2010. Meshless method for shallow water equations with free surface flow. Appl. Math. Comput. 217 (11), 5113–5124. Daoud, A.H., Rakha, K.A., Abul-azm, A.G., 2008. A two-dimensional finite volume hydrodynamic model for coastal areas: model development and validation. Ocean Eng. 35, 150–164. DHI, 2005. M21HD: Scientific Documentation, MIKE21 Flow Manual, Hydrodynamic module, Danish Hydraulic Institute, Horsholm, Denmark. Du, C., 2000. An element-free Galerkin method for simulation of stationary twodimensional shallow water flows in rivers. Comput. Method. Appl. M. 182 (1–2), 89–107. Gu, Y.T., 2005. Meshfree methods and their comparisons. Int. J. Comput. Methods 2 (4), 477–515. Gu, Y.T., Liu, G.R., 2005. Meshless methods coupled with other numerical methods. Tsighua Sci. Tec. 2, 8–15. Gutierrez, R.H., Laura, P.A.A., Rossi, R.E., 1994. The method of differential quadrature and its application to the approximate solution of ocean engineering problems. Ocean Eng. 21 (1), 57–66. Hanert, E., Roux, D.L., Legat, V., Deleersnijder, E., 2005. An efficient Eulerian finite element method for the shallow water equations. Ocean Model. 10 (1–2), 115– 136. Hashemi, M.R., Abedini, M.J., Malekzadeh, P., 2006. Numerical modeling of long waves in shallow water using Incremental Differential Quadrature Method. Ocean Eng. 33, 1749–1764. Hashemi, M.R., Abedini, M.J., Malekzadeh, P., 2007. A differential quadrature analysis of unsteady open channel flow. Appl. Math. Model. 31, 1594–1608. Hashemi, M.R., Abedini, M.J., Neill, S.P., Malekzadeh, P., 2008. Tidal and surge modeling using differential quadrature: A case study in the Bristol Channel. Coast. Eng. 55, 811–819. Homayoon, L., Abedini, M.J., Hashemi, M.R., 2013. RBF-DQ solution for shallow water equations. J. Waterw. Port C ASCE. 139 (1), 45–60. Hon, Y. C., Cheung, K. F., Mao, X. Z., Kansa E. J., 1999. Multiquadric solution for shallow water equations. J. Hydraul. Eng. ASCE 125 (5), 524–533. Kesserwani, G., Liang, Q., 2010. A discontinuous Galerkin algorithm for the twodimensional shallow water equations. Comput. Method. Appl. M. 199 (49–52), 3356–3368. Khoshfetrat, A., Abedini, M.J., 2012. A hybrid DQ/LMQRBF-DQ approach for numerical solution of Poisson-type and Burger’s equations in irregular domain. Appl. Math. Model. 36, 1885–1901. Rasulov, M., Aslan, Z., Pakdil, O., 2005. Finite differences method for shallow water equations in a class of discontinuous functions. Appl. Math. Comput. 160 (2), 343–353. Ren, Y., Wu, X., 2006. Boundary reduction technique and triangular differential quadrature domain decomposition method for polygonal region. Eng. Anal. Bound. Elem. 30, 435–440. Rodriguez-Paz, M., Bonet, J., 2005. A corrected smooth particle hydrodynamics formulation of the shallow-water equations. Comput. Struct. 83 (17–18), 1396– 1410. Shu, C., 2000. Differential Quadrature and its Application in Engineering. SpringerVerlag, London. Shu, C., Chew, Y.T., 1997. Fourier expansion-based differential quadrature and its application to Helmholtz eigenvalue problems. Commun. Numer Methods Eng. 13, 643–653.
10
A. Khoshfetrat, M.J. Abedini / Ocean Modelling 65 (2013) 1–10
Shu, C., Chew, Y.T., Liu, Y., 1998. Different interface approximations in multi-domain GDQ simulation of Czochralski bulk flows. Int. J. Numer. Method. H. 8 (4), 424– 444. Shu, C., Ding, H., Yeo, K.S., 2003. Local radial basis function-based differential quadrature method and its application to solve two-dimensional incompressible Navier–Stokes equations. Comput. Method. Appl. M. 192, 941–954. Shu, C., Richards, B.E., 1992. Application of generalized differential quadrature to solve two-dimensional incompressible Navier–Stokes equations. Int. J. Numer. Method. Fluid. 15, 791–798. Skiba, Y.N., Filatov, D.M., 2008. Conservative arbitrary order finite difference schemes for shallow-water flows. J. Comput. Appl. Math. 218 (2), 579–591. Toda, K., Ogata, Y., Yabe, T., 2009. Multi-dimensional conservative semi-Lagrangian method of characteristics CIP for the shallow water equations. J. Comput. Phys. 228 (13), 4917–4944. Wang, X., Liu, P.L.-F., 2011. An explicit finite difference model for simulating weakly nonlinear and weakly dispersive waves over slowly varying water depth. Coast. Eng. 58 (2), 173–183.
Wong, S.M., Hon, Y.C., Golberg, M.A., 2002. Compactly supported radial basis functions for shallow water equations. Appl. Math. Comput. 127 (1), 79–101. Wu, X.H., Kong, W.B., 2005. A highly accurate linearized method for free boundary problems. Comput. Math. Appl. 50, 1241–1250. Wu, X.H., Liu, S.T., 2005. Differential quadrature domain decomposition method for problems on a triangular domain. Numer. Meth. Part. D. E. 21, 574–585. Wu, X.H., Shen, Y., 2004. Differential quadrature domain decomposition method for a class of parabolic equations. Comput. Math. Appl. 48, 1819–1832. Wu, Z.M., 1992. Hermite–Birkhoff interpolation of scattered data by radial basis function. Approx. Theory Appl. 8, 1–10. Zienkiewicz, O.C., 1989. The finite element method. McGraw-Hill. Zong, Z., Lam, K.Y., Zhang, Y.Y., 2005. A multidomain differential quadrature approach to plane elastic problems with material discontinuity. Math. Comput. Model. 41, 539–553.
Contents lists available at SciVerse ScienceDirect
Ocean Modelling journal homepage: www.elsevier.com/locate/ocemod
Numerical modeling of long waves in shallow water using LRBF-DQ and hybrid DQ/LRBF-DQ A. Khoshfetrat 1, M.J. Abedini ⇑ Dept. of Civil and Environmental Engineering, Shiraz University, Shiraz, Iran
a r t i c l e
i n f o
Article history: Received 13 August 2012 Received in revised form 12 January 2013 Accepted 20 January 2013 Available online 15 February 2013 Keywords: Differential quadrature Local RBF-based differential quadrature Shallow water Oresund strait Hybrid methods
a b s t r a c t A computational algorithm based on coupling the differential quadrature (DQ) method and the local radial basis function-based differential quadrature (LRBF-DQ) method is used to solve the shallow water equations. The proposed method that is named ‘hybrid DQ/LRBF-DQ’ has the advantages of DQ such as simplicity and low computational cost and the advantages of LRBF-DQ, like mesh free and Kronecker delta function properties. Contrary to the DQ method, this method can be used in problems with irregular domain with less computational cost and less sensitivity to shape parameter(s) compared to the LRBF-DQ method. The proposed method is employed to simulate the complex and irregularly shaped waterway of the Oresund strait which is located between Sweden and Denmark. Implementation of the proposed method leads to good results for water depth. In comparison with the LRBF-DQ method, the results show that the proposed method is less sensitive to the shape parameter and the range of acceptable shape parameter in the proposed method is wider. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction The calculation of tidal elevations and currents in coastal regions constitute an important part of any coastal management study. As solving the governing partial differential equations (PDEs) analytically is impossible, an approximate numerical solution must thus be found using a suitable numerical method. The most popular methods are the finite difference (Rasulov et al., 2005; Skiba and Filatov, 2008; Wang and Liu, 2011), the finite element (Hanert et al., 2005; Cotter et al., 2009; Kesserwani and Liang, 2010), the finite volume (Castro et al., 2008; Daoud et al., 2008; Adamy et al., 2010) and the method of characteristics (Toda et al., 2009). The Differential Quadrature Method (DQM) proposed by Bellman and Casti (1971), Bellman et al. (1972) is an alternative numerical method for the numerical solutions of linear and nonlinear partial differential equations. Using concepts borrowed from polynomial approximation and linear vector space, Shu and Richards (1992) presented the polynomial-based differential quadrature (PDQ), in which the weighting coefficients of any order derivative are calculated by a simple algebraic formulation or a recurrence relationship without any restriction on the choice of grid distribution. Later, Shu and Chew (1997) also developed the Fourier-expansion-based differential quadrature (FDQ) where the ⇑ Corresponding author. Tel.: +98 711 6474604; fax: +98 711 6473161. 1
E-mail address: [email protected] (M.J. Abedini). Current address: Khorasgan Branch, Islamic Azad University, Isfahan, Iran.
1463-5003/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ocemod.2013.01.006
function is approximated by the Fourier series expansion. The details of PDQ and FDQ as well as their applications in different disciplines can be found in Shu (2000). Due to its efficiency and accuracy, DQM has been widely employed in many areas of industry and mathematics (Wu and Shen, 2004; Wu and Kong, 2005; Wu and Liu, 2005; Ren and Wu, 2006). Gutierrez et al. (1994) proposed the implementation of DQM in the field of ocean engineering. Hashemi et al. (2006, 2007, and 2008) used DQM as a rapid and accurate method for numerical simulation of one dimensional shallow water long waves. They showed the advantages for DQM compared with other conventional methods, namely, the accuracy, simplicity and minimal computational effort. For the implementation to become a practical alternative to existing integration methods, further research is necessary. As an example, the method should address the extension of the method to a tidal problem in an irregular domain. On the other hand, the use of mesh free (MFree) methods has been a subject of research interest in physical and engineering problems. In the MFree methods, the field function is constructed entirely from a set of arbitrarily distributed nodes (Belytschko et al., 1996; Gu, 2005). The nodal connectivity is imposed with the use of influence domains instead of elements being conventional in the finite element method (Zienkiewicz, 1989). These methods have been recently used to model shallow water equations (Du, 2000; Rodriguez-Paz and Bonet, 2005; Darbani et al., 2010; Hon et al., 1999; Wong et al., 2002). Homayoon et al. (2013) used RBF-DQ that is an MFree method, for modeling the shallow water equations. This method suffers
2
A. Khoshfetrat, M.J. Abedini / Ocean Modelling 65 (2013) 1–10
from two main drawbacks which are high computational cost and sensitivity to shape parameter. Shu et al. (2003) introduced Local RBF-DQ (LRBF-DQ) in which the two drawbacks of RBF-DQ are reduced. In this paper, LRBF-DQ method will be used to model shallow water equations in Oresund strait. It is often desirable and beneficial to combine two established numerical methods to exploit their advantages while avoiding their disadvantages (Gu and Liu, 2005). Khoshfetrat and Abedini (2012) introduced a new method, which is named hybrid DQ/LRBF-DQ, through coupling the conventional DQ and the LRBF-DQ to solve PDEs. For this purpose, the computational domain is divided into a few rectangular shapes and some irregular shapes. In such a domain decomposition process, a high percentage of the computational domain will be covered by regular shapes thus taking advantage of conventional mesh-based DQM eliminating the need to implement local RBF-DQ over the entire domain but only on a portion of it. By this method, we have the advantages of DQ like simplicity, high accuracy, and low computational cost and the advantages of LMQRBF-DQ like mesh free and Kronecker delta function properties. While irregular geometry has been the main impetus to use domain decomposition in this study, the multi-domain technique has been implemented in other studies for totally different purposes (Shu et al., 1998; Zong et al., 2005). Shu et al. (1998) used multi-domain technique to remove the effect of singularity and Zong et al. (2005) proposed domain decomposition to address material discontinuity in the computational domain. It is expected that the computational cost in the proposed methodology be less than the LRBF-DQ method. It will also be shown that the sensitivity to shape parameter in the proposed methodology is less than the LRBF-DQ method. We will use the hybrid DQ/LRBF-DQ method to model the shallow water equations in Oresund strait and the results will be compared with observed data, the output of MIKE21 model and that obtained by LRBF-DQ method. 2. Shallow water equations (SWEs) Shallow water free surface flows are usually described by the SWEs which can be derived from the depth averaged Navier– Stokes equations over the flow depth. Assuming the fluid density is constant and the pressure distribution is hydrostatic, the continuity equation can be written as:
@h @P @Q þ þ ¼0 @t @x @y
ð1Þ
u¼
1 h
Z
zs
Udz;
v¼
zf
1 h
Z
zs
Vdz
ð5Þ
zf
where U and V are the local velocity components and Zf and Zs are the bed and surface elevations. In Eq. (3) g is the acceleration of gravity; n is the Manning roughness coefficient; f = 2xSink is the Coriolis parameter; x is the angular velocity of the earth rotation; k is the latitude of the studied area; q is the water density; qair is the air density; Cw is the wind friction factor; Uw and Vw are the wind speed components in the x and y directions. The governing equations (1)–(3) are subject to appropriate initial and boundary conditions. Usually, the domain is limited by the bottom, the water surface, the land boundaries and the boundaries which could be purely imaginary that is named open boundaries. The boundary conditions on bottom and free surface have been considered in deriving the SWEs. Strictly speaking, the boundary condition for velocity on a land boundary is a no-slip condition: ~ ¼ 0, in which U ~ ¼ ðU; VÞ is the velocity vector. However, due to U the turbulence and the existence of a boundary layer, the velocity near a wall quickly becomes non-zero, and the no-slip condition is often replaced by the impermeability or free slip boundary condi~~ tion: U; n ¼ 0, where ~ n is the unit outward normal vector on the land boundary. On the open boundaries, additional information is required on the sea surface elevation as a function of time: Z s ¼ Z s ðtÞ, where Z⁄ is the specified sea surface elevation on the ~ ¼ ð0; 0Þ, which correopen boundary. The initial conditions are U sponds to the still water condition and Z s ¼ cte, that is the mean sea level in the computational domain. 3. Hybrid DQ/LRBF-DQ method 3.1. Description In the hybrid DQ/LRBF-DQ method, the domain of problem is decomposed into some regular and irregular subdomains. Consider a typical two-dimensional problem as shown in Fig. 1a. The problem domain consists of five subdomains X1 ; X2 ; X3 ; X4 ; X5 . If the point (xi, yj) is located inside region X1, then the conventional DQM will be used to approximate the derivative with respect to x, ux and the derivative with respect to y, uy at (xi, yj) via:
ux ðxi ; yj Þ ¼
NX X wxi;k uðxk ; yj Þ
ð6aÞ
k¼1
uy ðxi ; yj Þ ¼
NY X wyj;k uðxi ; yk Þ
ð6bÞ
k¼1
The momentum equation in x-direction can be expressed as: 2
!
@Z f @P P @h 2P @P PQ @h Q @P P @Q þ gh 2 þ þ þ þ gh @t h @x h2 @y h @y h @y @x h @x qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi gn2 P P2 þ Q 2 q þ fQ air C w U w U 2w þ V 2w ¼ 0 7=3 q h
ð2Þ
where NX, NY are, respectively, the number of mesh points in region X1 in the x and y directions, wxi,k, wyj,k are the DQ weighting coefficients in the x and y directions. In order to determine the above weighting coefficients, many kinds of test functions have been used including polynomials, the sine functions, the cosine functions and various orthogonal polynomials. The Lagrange interpolation function is widely used for this purpose, namely
For the y-direction, the momentum equation is:
lk ðxÞ ¼
2
@Z f @Q PQ @h Q @P P @Q Q @h 2Q @Q þ 2 þ þ þ ðgh 2 Þ þ þ gh @t @y h @y h @y @y h @y @y h h qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi gn2 Q P2 þ Q 2 q þ fP air C w V w U 2w þ V 2w ¼ 0 þ 7=3 q h
ð3Þ
where t is time; h is the flow depth; P and Q are the discharge per unit width in x and y directions that can be expressed as:
P ¼ uh;
Q ¼ vh
ð7Þ
Substituting Eq. (7) into Eq. (6a), we can get the following formulas to directly calculate the weighting coefficients of the firstand second-order derivatives, i.e.
wxi;j ¼
ð4Þ
In which, u and v are the depth averaged velocities in x and y directions which are defined respectively as:
NX Y ðx xi Þ x xi i¼1;i–k k
wxi;i ¼
NX Y 1 xi xk ; for j–i xj xi k¼1;k–i;j xj xk NX X
1 x xk k¼1;k–i i
ð8Þ
ð9Þ
A. Khoshfetrat, M.J. Abedini / Ocean Modelling 65 (2013) 1–10
(a)
3
(b)
Fig. 1. Domain decomposition for Hybrid DQ/RBF-DQ method.
The weighting coefficients in the y direction can be obtained in a similar way. If the point ðxl ; yl Þ is located in the irregular subdomains X2 ; X3 ; X4 ; X5 (containing the interface regions), local RBF-DQ is used to find the derivative with respect to x, ux and the derivative with respect to y, uy at ðxl ; yl Þ by:
ux ðxl ; yl Þ ¼
N X wxl;m uðxm ; ym Þ
ð10aÞ
m¼1
uy ðxl ; yl Þ ¼
N X wyl;m uðxm ; ym Þ
ð10bÞ
m¼1
where N is the number of supporting points in the neighborhood of ðxl ; yl Þ; wxl;m and wyl;m are the LRBF-DQ weighting coefficients and u(xm,ym) is the function value at all supporting points (xm,ym). The weighting coefficients in Eq. (10) are determined by using the multi-quadric (MQ) RBFs as test functions. MQ RBFs can be expressed as follows:
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi up ðXÞ ¼ uðjjX X p jjÞ ¼ uðrÞ ¼ r2 þ c2
ð11Þ
where X ¼ ðx; yÞ; r ¼ jjX X p jj is the Euclidean norm, and c is a free shape parameter. Substituting the set of MQ radial basis functions into Eqs. (10), one can obtain: N @ up ðxl ; yl Þ X ¼ wxl;m up ðxm ; ym Þ @x m¼1
ð12aÞ
N @ up ðxl ; yl Þ X wyl;m up ðxm ; ym Þ ¼ @y m¼1
ð12bÞ
For p ¼ 1; 2; . . . ; N. Equation systems (12a) and (12b) are used to compute the weighting coefficients wxl;m and wyl;m , respectively. For the case when the total number of points in the computational domain are used in Eqs. (10) and (12), then the derivatives is said to be approximated by the RBF-DQ method. Now, we can combine the DQ and LRBF-DQ equations used in hybrid DQ/LRBF-DQ method, into one system of equations:
ux ðxi ; yi Þ ¼
Nix X wxi;j uðxj ; yj Þ
ð13aÞ
j¼1
uy ðxi ; yi Þ ¼
Niy X wyi;j uðxk ; yk Þ
ð13bÞ
k¼1
where (xi,yi) is the point at which the derivatives are to be approximated. These are the formulas for derivative approximation in the proposed hybrid DQ/LRBF-DQ method. To determine the Nix, Niy and points (xj,yj), (xk,yk) in equations (13a, b), If the point (xi,yi) is located inside region X1, then the conventional DQM will be used to
compute the weighting coefficients in Eqs. (13a, b), so the Nix and Niy are, respectively, the number of rectangular mesh points (including the boundary points) in the x and y directions (Fig. 1b), and yj = yi, xk = xi in Eqs. (13a) and (13b), respectively. If the point (xi,yi) is located in the irregular subdomains X2,X3,X4,X5 (containing the interface regions), the Local RBF-DQ is used to find the weighting coefficients in Eqs. (13a) and (13b). Nix = Niy are the number of supporting points surrounding the point under consideration and uðxj ; yj Þ ¼ uðxk ; yk Þ is the function value at the supporting point. Supporting points may be located either in regular or irregular regions. Once again, spatial and temporal connectivity of state variable(s) between regular and irregular domains will be taken care of by points common to both regions. By using this method, we have the low cost DQM in a large part of domain and a mesh free Local RBF-DQ in the small regions near the boundaries for implementation of non-rectangular geometry which cannot be considered by the conventional DQM (Fig. 1b). 3.2. Advantages over traditional methods - The proposed hybrid approach and its variants (i.e., DQ and RBF-DQ) are categorized as higher order accurate due to reservation of a large number of terms in the Taylor series expansion. The exact order depends on the number of nodes considered to expand the derivative. - For a conventional triangular linear element extensively used in FEM and FVM, the state variable(s) would vary linearly between nodal points, while in RBF-DQ with multi-quadric basis function, the state variable(s) would vary in an exponential manner between nodal points. The major implication of variation in state variable(s) between nodal points concerns with rate of convergence. Multi-quadric RBFs have an exponential rate of convergence, while the conventional FEM and FVM implementation has linear rate of convergence. That is why it is argued that with much fewer nodes; the results obtained compared well with analytical and benchmark solutions of other numerical schemes. In other words, with a similar node distribution in RBF-DQ and FE, by similarly refining the nodes, the RBF-DQ results converge to the correct solution more rapidly than those of FEM or FVM. Furthermore, it might be interesting to know that one more advantage of RBF-based method as summarized by Wu (1992) is that it is space dimension independent in the sense that the convergence order is of O(hd + 1) where h is the density of the collocation points and d is the spatial dimension. - It is a well-known fact that traditional numerical techniques such as FDM, FEM, and FVM are strongly dependent on mesh properties. In light of this, mesh generation takes up very much time for complex geometries while resorting to these methods for obtaining numerical solution of continuum mechanics problems. Indeed, the mesh generation is considered to be the most expensive part of the simulation in these methods. Mesh-free methods (also called mesh-less methods) such as RBF-DQ have
4
A. Khoshfetrat, M.J. Abedini / Ocean Modelling 65 (2013) 1–10
generated considerable interest in recent years due to the need to overcome the high cost of mesh generation associated with the human labor in mesh-based methods. In mesh-free methods, such mesh generation is not required and instead, node generation is done which is much easier. In node generation, node connectivity and the area between the nodes is not important as opposed to in mesh generation. As an example, in FEM, if a triangular element becomes so narrow, poor results is obtained. However, in RBF-DQ, the characteristics of the area between the nodes are of no importance. As a result, node generation in RBF-DQ is very easier than mesh generation in FEM. Actually, the notion ‘‘mesh-free’’ is more suitable than ‘‘meshless’’, as in these methods, the nodes could have more degree of freedom compared to mesh-based methods. In short, the cost of mesh generation is a big issue in mesh-based methods. Generally speaking, in contrast to mesh-based methods where the grid resolution and/or choice of basis function are the primary determinant of accuracy, the shape parameter plays an important role in the accuracy of RBF-based methods in addition to the node distribution/density. The salient feature of our proposed hybrid numerical scheme is to diminish the impact of shape parameter on numerical accuracy. - In mesh-based methods such as FEM, the interpolating function will be defined over a typical element, while in mesh-free methods such as RBF-DQ; the interpolating function will be defined over a typical node as opposed to a typical element. The major implication of this change in focus from ‘‘element to node’’ relates to node refinement. In mesh-free methods, node refinement can be handled very easily, while in meshbased methods, re-meshing would be quite cumbersome for irregular domains if the numerical exercise requires doing so. - In FEM (or FVM), definition of interpolating function over a typical element is highly localized, while in RBF-DQ, its definition has a global nature. Indeed, the state variable has nonzero values over a typical element in FEM (or FVM) and zero value elsewhere. This won’t be the case in RBF- DQ while using infinitely smooth RBFs such as multi-quadric. The shape parameters are intended to control the impact of neighboring nodal points.
4. Numerical formulation of SWEs In this paper, the hybrid DQ/LRBF-DQ and LRBF-DQ methods are adapted to solve the shallow water equations (1)–(3) subject to the initial and boundary conditions as described in Section 2. Eqs. (1)– (3) can be approximated as: tþ1
hi
t
Nix X hi X wxij Ptj þ wyik Q tk ¼ 0 þ Dt j¼1 k¼1 Niy
ð14Þ
0 !2 1 Niy Nix Nix P tþ1 P ti @ t P ti AX 2P t X Pt Q t X t t i wxij hj þ t i wxij P tj i t 2i wyik hk þ ghi t Dt hi hi j¼1 ðhi Þ k¼1 j¼1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 t Niy Niy tX Nix X gn P ðP ti Þ þ ðQ ti Þ Q ti X P i t þ t wyik P tK þ ti wyik Q tk þ ghi wxij Z fj þ t 7=3 hi k¼1 hi k¼1 ðhi Þ j¼1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q t 2 2 ð15Þ fQ i þ air C w U tw ðU tw Þ þ ðV tw Þ ¼ 0
q
Nix Nix Nix Q itþ1 Q ti Pti Q ti X Qt X Pt X t wxij hj þ ti wxij Ptj þ it wxij Q tj t 2 Dt hi j¼1 hi j¼1 ðhi Þ j¼1 0 !2 1 Niy Niy Niy X Q ti AX 2Q t X t t t @ þ ghi wyik hk þ t i wyik Q tk þ ghi wyik Z fk t hi hi k¼1 k¼1 k¼1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi gn2 Q ti ðPti Þ2 þ ðQ ti Þ2 q t þ þ fPi air C w V tw ðU tw Þ2 þ ðV tw Þ2 ¼ 0 t 7=3 q ðhi Þ
ð16Þ where wxij and wyik are the hybrid DQ/LRBF-DQ or the LRBF-DQ weighting coefficients that are determined according to previous section. Nix and Niy are the number of supporting points in the LRBF-DQ method or the number of points in x or y direction in DQ method. The subscript and superscript of h, P and Q denote the location and time, in which the equations are approximated, respectively. It can be seen that the equations are discretized by explicit method in time, or the state variables in time step t + 1 are unknown.
Fig. 2. Location of Oresund strait (inside the circle).
A. Khoshfetrat, M.J. Abedini / Ocean Modelling 65 (2013) 1–10
~~ For satisfying the impermeability condition: U: n ¼ 0 on the land boundary, the flow velocity on this boundary needs to be adjusted so that the normal flow velocity is equal to zero. Here, we have the discharge per unit width instead of flow velocity. In each time step, the computed discharges on the land boundaries (P and Q) are updated as:
P ¼ Pn2y Qnx ny
ð17Þ
Q ¼ Qn2x Pnx ny
ð18Þ
where nx and ny are the x and y components of the normal vector on the land boundary, respectively. Eqs. (17) and (18) eliminate the normal component of the calculated discharge on the land boundary.
5
5. Solving SWEs in Oresund strait Oresund strait, located between Denmark and Sweden, connects Baltic Sea to Kattegat gulf (Fig. 2). In 1994, the construction of a fixed link between Copenhagen (Denmark) and Malmo (Sweden) as a combined tunnel, bridge and reclamation project commenced. Severe environmental constraints were enforced to ensure that the environment of the Baltic Sea remains unaffected by the link. To meet the environmental constraints and to monitor the construction work, a major monitoring program was set up. The resulting data were used to calibrate the DHI’s three-dimensional model, MIKE 3 and its two-dimensional model, MIKE 21. MIKE21 Hydrodynamic Software considers Alternating Direction Implicit (ADI) technique to integrate the equations for mass and momentum conservation in the space–time domain (DHI, 2005).
Fig. 3. The limits of modeling domain (inside the dashed rectangular).
6
A. Khoshfetrat, M.J. Abedini / Ocean Modelling 65 (2013) 1–10
Here, the data reported in MIKE 21 are used to evaluate the LRBFDQ and hybrid DQ/LRBF-DQ methods and the modeling results are compared with MIKE 21 results. The part of Oresund strait located in the 63900m 83700m dashed lines rectangular that is shown in Fig. 3 is modeled. The x and y axis are selected in a way that at open boundaries, the flow
direction is perpendicular to boundary of domain. The origin is located at Easting 337100 and Northing 6122900. The domain is discretized as shown in Fig. 4. The grid spacing in x and y directions is 900m. These grid points which are considered in MIKE 21 model, is used for modeling the Oresund strait by the LRBF-DQ and the hybrid DQ/LRBF-DQ methods.
Fig. 4. Grid points considered in modeling domain and the location of measuring stations.
Fig. 5. Comparison of LRBF-DQ and measured and MIKE21 simulated results at Drogden station.
A. Khoshfetrat, M.J. Abedini / Ocean Modelling 65 (2013) 1–10
In addition to the bathymetry of all grid points that is available, the time series of the following data has been measured in 1800 s time intervals in a 12 days period from 12/2/1993 to 12/ 13/1993:
7
- Water depth at two points of each open boundary (WL1, WL2, WL3 and WL4) in Fig. 4. These data are used to compute the water depth at all grid points in the open boundaries by linear interpolating method.
Fig. 6. Comparison of LRBF-DQ and measured and MIKE21 simulated results at NdrRoese station.
Fig. 7. Domain decomposition of Oresund strait for hybrid DQ/LRBF-DQ method.
8
A. Khoshfetrat, M.J. Abedini / Ocean Modelling 65 (2013) 1–10
- Water depth at Drogden and NdrRoese stations (Fig. 4) that are used for model validation. - Wind speed and direction at strait. The water elevation and discharge per unit width at t = 0 are considered 0.38m and zero, respectively (initial conditions). The LRBF-DQ method has been used to model SWEs in Oresund strait and the results are shown in Figs 5 and 6. In these figures, the computed water elevations at Drogden and NdrRoese stations during the 12 days measuring period are plotted. To evaluate the LRBF-DQ method in modeling SWEs, the measured water elevations and the results of MIKE 21 model are plotted in these figures as well. The time step of modeling is 100 s and the number of supporting points used to compute the shape functions is 200. The various values have been considered for the shape parameter and the best results are shown in Figs. 5 and 6. Implementation of the hybrid DQ/LRBF-DQ method involves the decomposition of computational domain into some regular and irregular subdomain. This decomposition for Oresund strait is shown in Fig. 7. It can be seen that each regular region contains four nodes in each of the x and y directions. It means that in each regular subdomain, the derivatives in four interior nodes are approximated by the low-priced DQ method. According to domain
decomposition shown in Fig. 10, in 34% of the nodes, the DQ method will be used for numerical modeling when using hybrid DQ/ LRBF-DQ method. This is why the experience on SWEs modeling in Oresund strait show that if the hybrid DQ/LRBF-DQ method is used, the computational time will be reduced to 84% of the computational time involved when the LRBF-DQ method is used. The computed time series of water depth by the hybrid DQ/LRBF-DQ method at Drogden and NdrRoese stations are plotted in Figs 8 and 9, respectively. The performance of hybrid method can be evaluated by comparing the numerical results, obtained from this method, with actual measurements and MIKE21 results. Here, in hybrid DQ/LRBF-DQ method, the modeling time interval is considered as 100 s and the number of supporting points is 200. The accuracy of results depends on shape parameter value selected. The various values have been considered for the shape parameter and the best results are plotted in Figs. 8 and 9. To assess and evaluate the performance accuracy of various approaches, the following error norm is considered:
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N N X X 1u Þ2 = u 2 Error ¼ t ðu u N i¼1 i¼1
Fig. 8. Comparison of hybrid DQ/LRBF-DQ and measured and MIKE21 simulated results at Drogden station.
Fig. 9. Comparison of hybrid DQ/LRBF-DQ and measured and MIKE21 simulated results at NdrRoese station.
ð19Þ
A. Khoshfetrat, M.J. Abedini / Ocean Modelling 65 (2013) 1–10
9
model and good agreement is achieved. The proposed method has the advantages cited in Section 3.2 compared to the MIKE21 model in which a traditional numerical method is used. The results show that the sensitivity of results to the shape parameter is reduced when using hybrid DQ/LRBF-DQ method. Additionally, the computational time is reduced in hybrid DQ/LRBF-DQ method. Results of this research show that the hybrid DQ/LRBF-DQ method is an attractive alternative for DQ and LRBF-DQ methods in ocean modeling problems that reduces their shortcomings.
References
Fig. 10. Variation of error with shape parameter (c) for various methods at Drogden station.
Fig. 11. Variation of error with shape parameter (c) for various methods at Drogden station.
where N is the number of data measured at stations, u is the simulated value of state variable obtained by a particular numerical is the measured data . method, and u To evaluate the influence of the value of shape parameter on accuracy of results obtained by the LRBF-DQ and the hybrid DQ/ LRBF-DQ methods, the variation of error with various values of shape parameter at Drogden and NdrRoese stations are plotted in Figs. 10 and 11, respectively. These figures show that the hybrid DQ/LRBF-DQ method is less sensitive to the shape parameter, for the reason that range of acceptable shape parameters in which the solution is convergent and steady, is larger and the error vs. c curve is less fluctuating.
6. Conclusions In this paper, the hybrid DQ/LRBF-DQ method is used in ocean modeling for the first time. This method prepares the opportunity of using the simple, accurate and cheap DQ method for problems with irregular domains, while the conventional DQ method can be used only in regular domains. Although the LRBF-DQ method is able to model problems with such irregular domains but it involves performing a large amount of computations. Another inconvenience in using LRBF-DQ is its sensitivity to the value selected for shape parameter. To illustrate the advantages of the proposed method, the SWEs are simulated in Oresund strait. The computed results by the LRBFDQ and the hybrid DQ/LRBF-DQ methods are compared with actual measurements and good agreement is indicated. The computed results are also compared with the results obtained by the MIKE21
Adamy, K., Bousquet, A., Faure, S., Laminie, J., Temam, R., 2010. A multilevel method for finite volume discretization of the two-dimensional nonlinear shallowwater equations. Ocean Model. 33(3–4), 235–256. Bellman, R., Casti, J., 1971. Differential quadrature and long term integration. J. Math. Anal. Appl. 34, 235–238. Bellman, R., Kashef, B.G., Casti, J., 1972. Differential quadrature: a technique for the rapid solution of nonlinear partial differential equations. J. Comput. Phys. 10, 40–52. Belytschko, T., Krongauz, Y., Organ, D., Fleming, M., Krysl, P., 1996. Meshless methods: an overview and recent developments. Comput. Methods Appl. Mech. Eng. 139 (1), 3–47. Castro, M.J., García-Rodríguez, J.A., González-Vida, J.M., Parés, C., 2008. Solving shallow-water systems in 2D domains using finite volume methods and multimedia SSE instructions. J. Comput. Appl. Math. 221 (1), 16–32. Cotter, C.J., Ham, D.A., Pain, C.C., 2009. A mixed discontinuous/continuous finite element pair for shallow-water ocean modelling. Ocean Model. 26 (1–2), 86–90. Darbani, M., Ouahsine, A., Villon, P., Naceur, H., Smaoui, H., 2010. Meshless method for shallow water equations with free surface flow. Appl. Math. Comput. 217 (11), 5113–5124. Daoud, A.H., Rakha, K.A., Abul-azm, A.G., 2008. A two-dimensional finite volume hydrodynamic model for coastal areas: model development and validation. Ocean Eng. 35, 150–164. DHI, 2005. M21HD: Scientific Documentation, MIKE21 Flow Manual, Hydrodynamic module, Danish Hydraulic Institute, Horsholm, Denmark. Du, C., 2000. An element-free Galerkin method for simulation of stationary twodimensional shallow water flows in rivers. Comput. Method. Appl. M. 182 (1–2), 89–107. Gu, Y.T., 2005. Meshfree methods and their comparisons. Int. J. Comput. Methods 2 (4), 477–515. Gu, Y.T., Liu, G.R., 2005. Meshless methods coupled with other numerical methods. Tsighua Sci. Tec. 2, 8–15. Gutierrez, R.H., Laura, P.A.A., Rossi, R.E., 1994. The method of differential quadrature and its application to the approximate solution of ocean engineering problems. Ocean Eng. 21 (1), 57–66. Hanert, E., Roux, D.L., Legat, V., Deleersnijder, E., 2005. An efficient Eulerian finite element method for the shallow water equations. Ocean Model. 10 (1–2), 115– 136. Hashemi, M.R., Abedini, M.J., Malekzadeh, P., 2006. Numerical modeling of long waves in shallow water using Incremental Differential Quadrature Method. Ocean Eng. 33, 1749–1764. Hashemi, M.R., Abedini, M.J., Malekzadeh, P., 2007. A differential quadrature analysis of unsteady open channel flow. Appl. Math. Model. 31, 1594–1608. Hashemi, M.R., Abedini, M.J., Neill, S.P., Malekzadeh, P., 2008. Tidal and surge modeling using differential quadrature: A case study in the Bristol Channel. Coast. Eng. 55, 811–819. Homayoon, L., Abedini, M.J., Hashemi, M.R., 2013. RBF-DQ solution for shallow water equations. J. Waterw. Port C ASCE. 139 (1), 45–60. Hon, Y. C., Cheung, K. F., Mao, X. Z., Kansa E. J., 1999. Multiquadric solution for shallow water equations. J. Hydraul. Eng. ASCE 125 (5), 524–533. Kesserwani, G., Liang, Q., 2010. A discontinuous Galerkin algorithm for the twodimensional shallow water equations. Comput. Method. Appl. M. 199 (49–52), 3356–3368. Khoshfetrat, A., Abedini, M.J., 2012. A hybrid DQ/LMQRBF-DQ approach for numerical solution of Poisson-type and Burger’s equations in irregular domain. Appl. Math. Model. 36, 1885–1901. Rasulov, M., Aslan, Z., Pakdil, O., 2005. Finite differences method for shallow water equations in a class of discontinuous functions. Appl. Math. Comput. 160 (2), 343–353. Ren, Y., Wu, X., 2006. Boundary reduction technique and triangular differential quadrature domain decomposition method for polygonal region. Eng. Anal. Bound. Elem. 30, 435–440. Rodriguez-Paz, M., Bonet, J., 2005. A corrected smooth particle hydrodynamics formulation of the shallow-water equations. Comput. Struct. 83 (17–18), 1396– 1410. Shu, C., 2000. Differential Quadrature and its Application in Engineering. SpringerVerlag, London. Shu, C., Chew, Y.T., 1997. Fourier expansion-based differential quadrature and its application to Helmholtz eigenvalue problems. Commun. Numer Methods Eng. 13, 643–653.
10
A. Khoshfetrat, M.J. Abedini / Ocean Modelling 65 (2013) 1–10
Shu, C., Chew, Y.T., Liu, Y., 1998. Different interface approximations in multi-domain GDQ simulation of Czochralski bulk flows. Int. J. Numer. Method. H. 8 (4), 424– 444. Shu, C., Ding, H., Yeo, K.S., 2003. Local radial basis function-based differential quadrature method and its application to solve two-dimensional incompressible Navier–Stokes equations. Comput. Method. Appl. M. 192, 941–954. Shu, C., Richards, B.E., 1992. Application of generalized differential quadrature to solve two-dimensional incompressible Navier–Stokes equations. Int. J. Numer. Method. Fluid. 15, 791–798. Skiba, Y.N., Filatov, D.M., 2008. Conservative arbitrary order finite difference schemes for shallow-water flows. J. Comput. Appl. Math. 218 (2), 579–591. Toda, K., Ogata, Y., Yabe, T., 2009. Multi-dimensional conservative semi-Lagrangian method of characteristics CIP for the shallow water equations. J. Comput. Phys. 228 (13), 4917–4944. Wang, X., Liu, P.L.-F., 2011. An explicit finite difference model for simulating weakly nonlinear and weakly dispersive waves over slowly varying water depth. Coast. Eng. 58 (2), 173–183.
Wong, S.M., Hon, Y.C., Golberg, M.A., 2002. Compactly supported radial basis functions for shallow water equations. Appl. Math. Comput. 127 (1), 79–101. Wu, X.H., Kong, W.B., 2005. A highly accurate linearized method for free boundary problems. Comput. Math. Appl. 50, 1241–1250. Wu, X.H., Liu, S.T., 2005. Differential quadrature domain decomposition method for problems on a triangular domain. Numer. Meth. Part. D. E. 21, 574–585. Wu, X.H., Shen, Y., 2004. Differential quadrature domain decomposition method for a class of parabolic equations. Comput. Math. Appl. 48, 1819–1832. Wu, Z.M., 1992. Hermite–Birkhoff interpolation of scattered data by radial basis function. Approx. Theory Appl. 8, 1–10. Zienkiewicz, O.C., 1989. The finite element method. McGraw-Hill. Zong, Z., Lam, K.Y., Zhang, Y.Y., 2005. A multidomain differential quadrature approach to plane elastic problems with material discontinuity. Math. Comput. Model. 41, 539–553.