Accurate and efficient hydrodynamic analysis of structures with sharp edges by the Extended Finite Element Method (XFEM): 2D studies

Achieving accurate numerical results of hydrodynamic loads based on the potential-flow theory is very challenging for structures with sharp edges, due to the singular behavior of the local-flow velocities. In this paper, we introduce the Extended Finite Element Method (XFEM) to solve fluid-structure interaction problems involving sharp edges on structures. Four different FEM solvers, including conventional linear and quadratic FEMs as well as their corresponding XFEM versions with local enrichment by singular basis functions at sharp edges, are implemented and compared. To demonstrate the accuracy and efficiency of the XFEMs, a thin flat plate in an infinite fluid domain and a forced heaving rectangle at the free surface, both in two dimensions, will be studied. For the flat plate, the mesh convergence studies are carried out for both the velocity potential in the fluid domain and the added mass, and the XFEMs show apparent advantages thanks to their local enhancement at the sharp edges. Three different enrichment strategies are also compared, and suggestions will be made for the practical implementation of the XFEM. For the forced heaving rectangle, the linear and 2nd order mean wave loads are studied. Our results confirm the previous conclusion in the literature that it is not difficult for a conventional numerical model to obtain convergent results for added mass and damping coefficients. However, when the 2nd order mean wave loads requiring the computation of velocity components are calculated via direct pressure integration, it takes a tremendously large number of elements for the conventional FEMs to get convergent results. On the contrary, the numerical results of XFEMs converge rapidly even with very coarse meshes, especially for the quadratic XFEM.


Introduction
Numerical analysis is playing an increasingly important role in marine hydrodynamics. Computational Fluid Dynamic (CFD) models based on the Navier-Storkes (NS) equations with proper turbulence modeling are the most comprehensive ones for this purpose. They are applicable in more applications than a potential-flow model, in particular when viscous flow separation and wave breaking become relevant and important. The computational costs, however, are normally too high to afford, which is regarded as one of the bottlenecks of CFD models, if they are heavily involved in the design of marine structures. Due to large-volume nature of most of the marine structures, the inertial effect is predominant whereas viscosity effect plays a secondary role. Therefore, the potential-flow theory is often applied together with empirical corrections for viscous effects.
For the potential-flow problems, Boundary Element Method (BEM) is the most commonly used numerical method in marine hydrodynamics, as it can reduce the dimension of the problem by one and only the boundaries of the fluid domain need to be discretized. Even though the number of unknowns is reduced in BEM compared with a volume method, it is still challenging for a conventional BEM to solve the resulting linear system with a large number of unknowns, because the matrix is dense. O(N 2 ) memory is required by the conventional BEMs, and O(N 2 ) and O(N 3 ) operations are required for iterative solvers and direct solvers, respectively. Here N denotes the number of total unknowns on the boundary surfaces.
Although BEM is a very popular numerical method in potential-flow hydrodynamic analyses, field solvers are also widely used. Wu and Eatock Taylor (1994) is among the first to use FEM to investigate 2D nonlinear free-surface flow problems in the time domain. Wu and Eatock Taylor (1995) studied the fully-nonlinear wave-making problem by both FEM and BEM, and suggested that FEM is more efficient than BEM in terms of both CPU time and computer memory. Ma et al. (2010a,b) used a FEM to simulate the interaction between 3D fixed bodies and steep waves. On the other hand, high-order volume methods have gained great interest. Bingham and Zhang (2007) and Engsig-Karup et al. (2009) developed 2D and 3D highorder Finite Difference Methods (FDMs) to study fullynonlinear water wave problems in potential flows. Shao and Faltinsen (2012) and Shao and Faltinsen (2014) proposed high-order Harmonic Polynomial Cell (HPC) methods in 2D and 3D respectively to study water waves and their interaction with structures. Some recent extensions were made to utilize immersed boundary strategies and overset meshes to achieve better accuracy and stability (e.g., see Hanssen et al., 2018;Tong et al., 2019Tong et al., , 2021Law et al., 2020;Liang et al., 2020). Compared to the BEMs, field solvers deal with sparse matrices, and the computational costs are roughly linearly dependent on the number of unknowns.
Ordinary boundary-element and volume methods, e.g. BEM, FEM, FDM and HPC methods, are based on local approximations using smooth functions. Thus, very fine meshes have to be applied at areas where the fluid solution tends to be singular. Sharp edges are widely present in typical offshore structures. Examples are pontoons of semi-submersibles and tension leg platforms (e.g., see Chen et al., 1995;Zhou and Wu, 2015), damping plates on offshore platforms (e.g., see Tao et al., 2007;Shao et al., 2016Shao et al., , 2019 and offshore floating wind turbine structures (e.g., see Xu et al., 2019), as well as the bilge keels on the ships. Besides, the analytical methods, such as the multiterm Galerkin method (e.g., see Li et al., 2019;Porter, 1995), have also been used to include the local singularities. From industrial application point of view, it is essential to be able to obtain accurate numerical results with affordable computational efforts. However, this is not always possible, even for the 2nd order mean wave loads.
The calculation of 2nd order mean wave loads involves quadratic terms of the 1st order quantities, which pose great challenges at the sharp edges where the fluid velocities tend to be infinite. Taylor and Teng (1993) investigated the effect of corners on diffraction/radiation wave loads and wave-drift damping, and revealed that the most important hydrodynamic loads and the amplitudes of body motion do not change significantly while the radius of the corner approaches zero. For a floating truncated vertical cylinder free to surge and heave, Zhao and Faltinsen (1989) found it is difficult to obtain convergent 2nd order mean wave forces via the direct pressure integration. In their work, a method based on momentum and energy relationship was shown to be more robust and efficient. By applying the variants of Stokes's theorem, Dai et al. (2005) and Chen (2007) developed a 'middle-field formulation', which transforms the body-surface integral to a control surface at a distance from the body. Similar strategy was also applied by Liang and Chen (2017) where a multidomain approach was developed. The middle-field formulation can be used to calculate drift forces and moments in all 6 degrees of freedom. The floating truncated vertical cylinder studied by Zhao and Faltinsen (1989) was revisited in Shao (2019) and four different methods were used to calculate the vertical mean wave force, including a momentum formulation implemented in a time-domain higherorder BEM (Shao and Faltinsen, 2013), a semi-analytical solution (Mavrakos, 1988), the middle-field method in Hy-droStar, as well as the near-field method in HydroStar.
The first three methods matched very well with each other, confirming the accuracy of the earlier results by Zhao and Faltinsen (1989) based on momentum and energy relations. However, the results determined by the direct pressure integration were quite different in the heave resonance regime. As elucidated in Shao (2019), the results by the direct pressure integration are not convergent, despite very fine meshes have been used. Yang et al. (2020) used five different methods to investigate nonlinear radiation forces of bodies with sharp or rounded edges in the time domain. The first four methods are all near-field methods, and the fifth one based on momentum conservation. They found that the singularity at the sharp edge plays significant roles on numerical computation of hydrodynamic forces in all near-field methods, while it has much less influences on results based on momentum conservation. Using an approach based on a control surface, Cong et al. (2020) rewrote the integration of velocity square terms on body surface into the sum of two other integrals, one on a control surface enclosing the structure and the other on the free surface between the structure and the control surface. Encouraging results were obtained for double-frequency wave-radiation forces on an oscillating truncated vertical cylinder.
This paper aims to introduce, verify and demonstrate the XFEM as an accurate and efficient tool to calculate the linear and 2nd order wave loads on structures with sharp edges, without having to use a control surface. The XFEM has a powerful framework, which allows for adding the knowledge of the local solutions, normally known as a priori, to the finite-element approximation space at specific nodes. The solution enrichment at those nodes does not require any modification to the meshes. The idea of XFEM was originally used by Belytschko and Black (1999) to solve the problem of elastic crack growth, and one year later, Daux et al. (2000) formally named this approach as XFEM. The XFEM can be seen as an extension of the standard FEM based on the conception of Partition of Unity (PU) (Babuška and Melenk, 1997), and thus it maintains all advantages of the standard FEM. Earlier concepts of PU dates back to 1994, when it was first used to solve the socalled roughness coefficient elliptic boundary value problem by Babuška et al. (1994) with the name of special finite element method, namely the Generalized Finite Element Method (GFEM). Based on the ideas in Babuška et al. (1994), the GFEM was further elaborated by Melenk (1995), Melenk andBabuška (1996), ? andMelenk andBabuska (1997) with the name of partition of unity method (PUM) or partition of unity finite element method (PUFEM). The GFEM was developed in Strouboulis et al. (2000a) and Strouboulis et al. (2000b) with the name of GFEM. In the early days, both XFEM and GEFM were developed independently even their basic idea is similar. A feature to distinguish the XFEM and the GFEM in early work is that only local parts of the domain are enriched by XFEM, but GFEM enriches the whole domain globally. However, Fries and Belytschko (2010) argued that the XFEM and the GFEM are almost identical numerical methods. The XFEM represents the singular properties by adding singular basis function or any analytical recog-nition of the solution to local approximation space, and it has been a tremendous success in dealing with singular or discontinuous problems, no matter how strong the discontinuity is (see, e.g. Sukumar et al., 2001;Moes et al., 2002;Sukumar et al., 2000;Fries and Belytschko, 2010). Besides, XFEM has also been introduced to CFD to model two-phase flows (Fries, 2010).
In the present work, as verification and demonstration, flow around an infinite-thin flat plate and a heaving rectangle on the free surface will be studied via four different FEMs, namely the linear FEM, linear XFEM, quadratic FEM and quadratic XFEM. Convergence studies will be presented to illustrate accuracy and efficiency of the XFEMs. Our results indicate that the singularities at sharp edges do not have a strong influence on the calculating of added mass and damping, confirming the conclusion from an earlier study by Taylor and Teng (1993). However, if the near-field method is used, it is extremely challenging for conventional FEMs to achieve convergent 2nd order vertical mean forces for the heaving rectangle with affordable computational time on a normal PC. On the contrary, the XFEMs with local enrichment, using cornerflow solutions (Newman, 2017) at sharp edges, can achieve convergent results with much coarser meshes. Three different local enrichment strategies of XFEM will also be compared and suggestions will be made for practical implementation.
The rest of the paper will be organized as follows. In Sect. 2, the formulation of the boundary-value problem and corner-flow solutions are presented. In Sect. 3, the basic idea of conventional FEM and the XFEM are introduced via a mixed boundary-value problem in 2D. Besides, three enrichment strategies for the XFEM are presented and compared. In Sect. 4, as the first verification case, the velocity potential in fluid domain and the added mass of an infinitely-thin flat plate are studied and compared with the analytical solution. The second verification concerns a heaving rectangle on a free surface, solved in the frequency domain. In Sect. 5, some conclusions are drawn.

Mathematical formulation
2.1. Governing equation and linearized boundary condition A 2D coordinate system Oxy is defined with the Ox axis coinciding with the undisturbed free surface and Oy axis orienting positively upward, as illustrated in Fig. 1. The fluid domain Ω is enclosed by the body surface S b , free surface S f , bottom surface S d , and vertical control surfaces S m at a distance from the body.
It is assumed that the fluid is inviscid and incompressible, and the flow is irrotational so that a velocity potential φ exists. In this study, we only consider 2D flows, and thus the governing equation in the fluid domain Ω is written as where φ denotes velocity potential. Here only radiation problem is considered, and thus the impenetrable condi-tion on the body surface is written as: where v is the velocity of the body and n is the vector normal to the body surface pointing out of the fluid domain. Besides, the combined linearized free-surface condition is written as The bottom condition is ∂φ

Linearized frequency-domain analysis
Assuming that the problem is time-harmonic and a steady state is reached. Therefore, velocity potential can be separated into a spatial part and temporal part as follows: where ω denotes the angular frequency of oscillation, and i = √ −1. The motion of body in j-th mode can be defined as: where η ja denote the amplitude of body in j-th mode, and j = 1, 2, and 3 correspond to sway, heave, and roll motions, respectively. Accordingly, the governing equation and boundary-value problem (BVP) with respect to the complex velocity potential ϕ(x, y) can be written as: Here n j represent the component of the normal vector in the direction of the motion of body in j-th mode. The dispersion relation in finite water depth is k tanh kh = ω 2 /g, where k is wavenumber. Thus, the free-surface condition in Eq. (7) can be rewritten as −(k tanh kh) · ϕ + ∂ϕ ∂y = 0 at y = 0.
Besides, radiation condition requiring radiated waves propagating outwards can be expressed as: If the horizontal distance between rectangle and matching boundary is large enough, the radiation condition can be satisfied at the matching boundaries S m :

Corner-flow solution
In order to demonstrate the singular characteristics of the corner flow by potential-flow theory, the flow past a sharp corner with an exterior angle β and the corresponding interior angle of γ = 2π − β as shown in Fig. 2 is considered. If the considered semi-infinite wedge is fixed, the corner-flow solution can, according to Newman (2017), be defined in the polar coordinate system Orθ as where A j is a constant and j is an non-negative integer number. It is obvious that the velocity determined by Eq. (11) is singular at the tip of the semi-infinite wedge when j ≥ 1 and γ < π. If we define Eq. (11) can be rewritten as For a general 2D radiation-diffraction problem, the local scatter velocity potential (incident wave potential excluded) close to a sharp edge can be expressed as Here the first term ϕ 0 is a smooth velocity potential satisfying the non-trivial Neumann-boundary conditions for the scatter velocity potential where ϕ i is the incident wave potential, v is the rigidbody velocity at the edge and n is the normal vector on the body surface. For a radiation problem, an example of ϕ 0 has been given by Liang et al. (2015) as ϕ 0 = u x + v y. u and v are the horizontal and vertical velocities at the corner due to rigid-body motions, respectively. The second term of Eq. (14) represents the corner-flow solutions derived from zero Neumann-boundary condition at the sharp edges. Since the first term is non-singular, it can be wellapproximated by ordinary shape functions. The singular terms in the second term (with j ≥ 1) are included in the local enrichment of XFEM to capture the local singular behavior at the edges.

Numerical method
Since the XFEM is an extension of the conventional FEM by including singular basis function in the shape function, we will in this section start with very brief introduction of the conventional FEM, followed by more details of XFEM as well as different local enrichment strategies for XFEM. General description of conventional FEMs can be found in many textbooks (see, e.g. Zienkiewicz et al., 2005;Hughes, 2012;Reddy, 2019).

Finite Element Method
In a FEM formulation for a potential-flow problem, the fluid domain is discretized into elements, also called finite elements, and the velocity potential in each element can be approximated as Here N j (x, y) is shape function, n p is the number of the nodes in the whole fluid domain and ϕ j denotes the nodal value of the velocity potential at node j. Application of the Galerkin method leads to For a mixed Dirichlet-Neumann BVP, as we will study in Sect. 4.1 for the flat plate in infinite fluid, ϕ = f p on S D and ∂ϕ/∂n = f n on S N are known from the the boundary conditions, respectively. In this case, Eq. (18) will be represented by a linear system as where Here the superscript T represents the transpose of a matrix or a vector. The elements in matrix K and vector B are defined respectively as For a mixed Neumann-Robin BVP, as we will study in Sect. 4.2 for the linear frequency-domain solution of a heaving rectangle at the free surface, the weak form can be more specifically written as: Here the mean free surface S f and control surface S m are Robin boundaries, where the boundary conditions are defined in Eqs. (8) and (9), respectively. The Neumann boundary condition on S b has been defined in Eq. (7).

Extended Finite Element Method (XFEM)
XFEM was developed based on the concept of partition of unity (PU), and the so-called PU means a set of non-zero function N i (x, y) in the partition of unity domain satisfying the following condition: For any function in the PU domain, the following relationship holds: Undoubtedly, Eq. (25) is also satisfied when ψ(x, y) is a constant. Obviously, standard shape functions, for instance those shown in Eqs. (48) and (49), are PU functions.
After introducing the conventional FEM and the conception of PU, the enrichment function and extra degrees of freedom (DOFs) at the selected nodes will be presented. For simplicity and without losing generality, we denote I as the set of all nodes in the fluid domain and J as the subset of nodes which will be enriched. Thus the trial solution in the fluid domain with only one enrichment function on each point j ∈ J can be written as where Ψ j represent the additional DOF at the enriched node j. N j (x, y) is the standard finite-element shape function, ψ(x, y) denotes the enrichment function representing special knowledge, e.g. logarithmic singularity, of the fluid solution. The products N j (x, y)ψ(x, y) may be considered as local enrichment function, as their supports coincide with those of conventional finite-element shape functions, leading to sparsity in the discrete equation (Fries and Belytschko, 2010). It can be understood from Eq. (26) that the values on nodes j ∈ J differ from ϕ j , which is an unfavorable property. To ensure that nodal values are always ϕ j at the enriched nodes j ∈ J , the enrichment function can be shifted and Eq. (26) can be rewritten as (see, e.g. Fries and Belytschko, 2010;Daux et al., 2000), As a result of the shifting, the enrichment represented by the 2nd summation on the right-hand side of Eq. (27) vanishes at the nodes j ∈ J , and thus recover the Kronecker-δ property of standard finite-element approximations. Unless otherwise redefined, all the enrichment functions that will be used in this paper are the shifted enrichment functions. More generally, if more than one enrichment function are introduced at each node j ∈ J , Eq. (27) can be extended as 5 Here ψ l (x, y) is the l-th enrichment function, ψ l (x j , y j ) denotes the value of ψ l (x, y) at j-th node, ψ l (x, y)−ψ l (x j , y j ) denotes the shifted enrichment function with a shifted value of ψ l (x j , y j ). For brevity, we use a matrix form to express Eq. (28) and rewrite it as where are the elements in N std and N enr . n enr denotes the number of enrichment functions. The dimension of N std is 1 × n p . If there are n enr p nodes enriched in whole domain, the dimension of N enr is 1 × (n enr p · n enr ). Substituting Eq. (29) into Eq. (18), we obtain the following expression: Here N D std denotes the shape function which lies on Dirichlet boundary, and Φ D represents the velocity potential of the nodes which are located on the Dirichlet boundary. We must emphasize that there are not any enrichment nodes on the Dirichlet boundary. For a mixed Dirichlet-Neumann BVP, in the same manner as Eq. (19), the linear system comes from Eq. (30) can be written as: The coefficient matrix of K can be divided into four parts as follows: where the elements in K ϕϕ , K ϕψ , K ψϕ and K ψψ are K ϕϕ comes from conventional standard finite elements, which is only relevant to standard shape function. K ϕψ , K ψϕ and K ψψ are related to the enrichment. X is a (n p + n enr p · n enr ) × 1 vector. The right-hand-side vector in Eq. (31) is where the elements in B ϕ and B ψ are B ϕ is related to the standard FEM and B ψ is related to the local enrichment. To be clear, Eq. (33) is equivalent to Eq. (21), and Eq. (38) is equivalent to Eq. (22) in the previous section. For a mixed Neumann-Robin BVP, we use corner-flow solution as the enrichment solution and Eq. (28) to construct the local approximation, and obtain a final equation system similar to Eq. (23) as: (40)

Enrichment strategies
In the previous subsection, the XFEM has been introduced through mixed BVPs. The key point of XFEM is the local enrichment, and we will in this subsection discuss three different enrichment strategies in detail. Unless otherwise mentioned in present work, a singular point is defined as a point where the fluid velocity is infinite, and an element is called singular element if it contains at least one singular point. A singular patch is a patch of multiple elements, among which at least one is a singular element.

Point enrichment
In the point enrichment approach, as depicted in Fig. 3, singular solutions are introduced to enrich the local approximation only on the singular points, and the end points of the blue line are the singular points. In this way, the additional number of unknowns due to enrichment is only dependent on the number of singular points and the number of singular terms introduced at each singular point, and thus is independent on the meshes. Therefore, this enrichment only influences the singular elements which contain the singular points. The influence domain of a enriched point depends on the mesh size. Consequently, the enhancement of the solution accuracy may not increase as the mesh are refined, which will be discussed later in Sect. 4.1.

Patch enrichment
Compared with point enrichment, the patch enrichment method will introduce enrichment to all points on the singular patch, which are represented by the filled circles in Fig. 4. The patch enrichment includes the first layer of neighboring points surrounding the singular points. Similar to the point enrichment method, the enrichment domain of patch enrichment depends on the mesh size, and the additional number of unknowns do not increase even if the meshes are refined. Similar to the point enrichment, patch enrichment experiences low convergent rate with the refinement of the local meshes.

Radius enrichment
Different from the point enrichment and patch enrichment, the radius enrichment method will enriches the solution at all points within a circle with predefined radius R enri . Here R enri must be a positive value and it is independent of the mesh size. As demonstrated in Fig. 5, the center point of the enrichment area is the singular point. The value of R enri may be taken as 1/10 of the space dimension, as it was suggested by (Laborde et al., 2005). In the present work, we normally take R enri = 0.2 as we are considering 2D problems. More detail about how to choose the enrichment radius will be discussed in the numerical example of heaving rectangle at free surface. The drawback of this enrichment method is that the additional number of unknowns will increase with the mesh refinements.

Integrals on the singular elements
In this subsection, the integration on singular elements is discussed. For illustrating the element integral strategy explicitly, a singular function stems from corner-flow solution in Sect. 2.3 is utilized as the enrichment function. In Eq. (13), the most singular term is the first term, i.e. the term with j = 1 ψ(x, y) = r m1 cos(m 1 θ).
For demonstration purpose, we will take this term as an example of the enrichment function, and discuss numerical integration of singular terms on the elements. In practice, more terms in Eq. (13) can be included as enrichment functions, using the similar procedure that will be described in the rest of this section. The trial solutions in Eqs. (26) and (27) involve the evaluation of the following enrichment shape function Here (x, y) is the location in physical space, which can be obtained from isoparametric element illustrated in Fig. 22. Derivatives of the enrichment shape function with respect to x and y are expressed by Substituting Eq. (43) into Eq. (36), the diagonal entry of enriched element stiffness matrix K ψψ ii can be written as where Ω e denotes the surface of elements. Point i is one of their nodes. Apparently, if the interior angle γ < π, the x-and y-derivatives of r m1 cos(m 1 θ) are singular, with a singularity of r m1−1 as r → 0. Thus the square terms in Eq. (44) are r 2m1−2 singularities. It is challenging but important to accurately calculate this singular integration.
In this paper, the so-called DECUHR adaptive quadrature algorithm (Espelid and Genz, 1994) is employed to overcome the difficulties in numerical integration of Eq. (44). The DECUHR algorithm combines an adaptive subdivision strategy with an extrapolation of the error expansion, where a non-uniform subdivision of the element close to a singular point is employed. More details of the DECUHR algorithm can be found in Espelid and Genz (1994), and the application of this algorithm in GFEM to deal with singular integrals can be found in Strouboulis et al. (2000a). An open-source FORTRAN code of this DECUHR algorithm from Alan Genz Software website of Washington State University, which can deal with problems at the dimension of 2∼15, has been applied in this study. It is not an option in the code to handle 1D singular integrals. In present work, an adaptive Gaussian quadrature algorithm is applied to accurately calculate the 1D singular integrals. It consists of the following steps: Step 1: Setting a fixed tolerance, using T to represent the result, letting T = 0, using T 3 represent the temporary variate and T 3 = 0.
Step 2: Using Gaussian integral to obtain the numerical integration result in whole element and written as T 1 , letting T = T 1 .
Step 3: Dividing the element into two uniform element, and integrating in those two subdivision element, expressing as T 21 and T 22 respectively, assuming T 21 is the numerical result in singular element, letting T 2 = T 21 + T 22 and T 3 = T 3 + T 22 .
Step 4: Calculate error = |T 2 −T 1 |. If error > tolerance, we denote T = T 3 + T 21 , divide the sub-element which contains singular point into two, and go to Step 2. Otherwise, if error ≤ tolerance, we output T as the final result.

Numerical studies
For verification purposes, an uniform flow around a flat plate of vanishing thickness in 2D, and the added mass of the same plate are considered. Then, the heaving rectangle on a free surface is studied via linear FEM, linear XFEM, quadratic FEM and quadratic XFEM.
4.1. Uniform flow around a flat plate and added mass of a flat plate The analytical solution of the complex potential for uniform flow around a 2D thin flat plate in an infinite domain can be found in the textbook of Newman (2017) and the Appendix B, where we also show that a modification of the sign in the original formula is needed for the flow variable on the right-half plane.
To model the flat plate in infinite fluid, we have to use a truncated fluid domain in our numerical method. See a sketch of the truncated domain in Fig. 6. Based on the analytical solution, Dirichlet boundary conditions are specified at the truncated boundaries surrounding the fluid domain, and Neumann boundary condition on the upper and lower surfaces of the thin plate. The mathematical formulation of the mixed BVP has been described in Sect. 2 and the conventional FEM and XFEM explained in Sect. 3. Even though the flat plate has zero thickness, the velocity potentials are different on the two sides of the plate. Thus a double-node technique is used on the plate except at the two endpoints of the flat plate, where the velocity potential must be continuous. The double-node technique allows for two velocity potential values at the same location. See an illustration in Fig. 6, where the open circles and crosses represent two different nodes respectively.
To solve the mixed Dirichlet-Neumann BVP numerically, we have implemented four different FEM solvers, including linear FEM, linear XFEM, quadratic FEM and  quadratic XFEM. In the linear and quadratic XFEMs, we have used the analytical solution as the enrichment function at the enrichment nodes close to the singular points, in this case the two ends of the plate. Figs. 3, 4 and 5 illustrate the point, patch and radius enrichment strategies, respectively. In order to compare the global accuracy of different methods, the L 2 errors of the velocity potential on all grid points will be presented as a function of mesh size ∆h = ∆x = ∆y. The L 2 error is defined as where φ num i denotes numerical solution on the ith node, and φ ana i represents the corresponding analytical solution. N denotes the total number of nodes.
In principle, the added mass of the plate should be calculated by considering oscillatory motions of the plate. However, for the special case of a structure in unbounded fluid or practically sufficient away from any other boundaries, the added mass can also be obtained based on the flow solution for fixed structure in a uniform flow. According to Section 4.10 in Newman (2017), the velocitypotential solution around a fixed structure in a free stream along y-axis can be expressed as φ fix = U y − φ move , where φ move = U φ 2 is the velocity potential due to the same structure moving along y-axis with velocity U . φ 2 is the velocity potential induced by the structure at a unit velocity, i.e. U = 1, along y-axis. Therefore, φ 2 can be immediately obtained as φ 2 = (U y − φ fix ) / U , and thus the added mass can be calculated according to Eq.(114) in Newman (2017), which only needs the flow solution of φ 2 . Note that the above discussions are still valid if the velocity U is not a constant in time, i.e. U = U (t).
The L 2 errors for velocity potential and the relatively errors for added mass of different method via different enrichment strategies are presented in Fig. 7, 8 and 9, respectively.
As shown in Fig. 7(a) and 7(b), when the point enrichment as described in Sect. 3.3.1 is applied at the two end-points of the plate, the errors are greatly reduced, indicating that the local enrichment in XFEMs is very effective in reducing both the local and global errors, which is   Figure 10: The error of velocity potential versus the non-dimensional enrichment radius R enri /a for linear XFEM. R enri represents the enrichment radius, a is half breadth of the plate. 64 × 64 uniform meshes have been used.
expected. However, it is also surprising that all FEMs, including the quadratic FEMs, showed convergent rates close to 1.0 for the velocity potential on fluid points. The k values in the figures are the fitted convergence rate based on five different mesh densities. Similarly, as seen in Fig. 7(b), lower than expected mesh-convergence rates are observed for the added mass of the flat plate. The influence area of the enrichment functions is smaller for a locally finer mesh close to the singular point, due to the fact that the interpolations by the finite-element shape functions in Eqs. (27) and (28) are piecewise. The shape function N j at the point j is always zero over the elements, which do not own point j as one of their element nodes. At the non-enriched points, sufficiently close to the singular point but belong to none of the singular elements, the velocity potential also changes dramatically, and thus the applied smooth shape functions have difficulties to accurately capture the strong local singular solutions. Actually, according to Zienkiewicz et al. (2005), the singularity affects not only the local area, but also at a distance surrounding the singularity, and consequently the convergence rate for the ordinary FEMs. The affected convergence rate follows O(N dof ) (−min[λ,p]/2) , where λ is a number associated with the intensity of the singularity, p represents power exponent of the basis function, N dof the number of freedom. More details about the effect can be found in Zienkiewicz et al. (2005, pp. 458-459).
The results of convergence studies are presented in Fig. 8 for the patch enrichment. For both the solutions of velocity potential at grid nodes and added mass, only marginal increases of the convergence rate of linear XFEM are seen. However, the improvement in the results of quadratic XFEM is notable for both velocity potential and added mass. Theoretically, we expect the convergence rate of a quadratic method be equal or greater than 2. Even though the overall accuracy of linear and quadratic XFEM has been greatly improved with the adoption of patch enrichment instead of the point enrichment, it is still below our expectation, in particular for the quadratic   XFEM. The reason is as follows: similar to the point enrichment strategy, the enrichment area of the patch enrichment strategy is also mesh-dependent.
To eliminate the mesh-dependency of the local enrichment, the radius enrichment method, as illustrated in Fig. 5, appears to be a good choice. In this method, the enrichment area is a predefined constant and independent on mesh sizes. As demonstrated by the results of mesh-convergence study in Fig. 9, the superiority of XFEMs, in particular the quadratic XFEM, is remarkable, when the radius enrichment is applied. We have used an constant enrichment radius of R enri = 0.2, as suggested by Laborde et al. (2005) for 2D problems, at both ends of the plate. Compared to the conventional linear FEM, convergence rate of linear XFEM for the velocity potential advanced notably from k = 0.89 to k = 1.38, and from k = 1.02 to k = 1.43 for added mass. The convergence rate of quadratic XFEM improved exceedingly from k = 0.99 to k = 3.44 for velocity potential, and for added mass from k = 1.11 to k = 1.79.
For the present case, since we are using the analytical solution as the enrichment function at the singular points, the accuracy of the XFEM solutions will further improve if a larger enrichment radius is applied. This is illustrated in Fig. 10, where we have presented the L 2 errors of the velocity potential as function of R enri /a, the ratio between enrichment radius and length of the plate.
For both linear and quadratic XFEMs, it is apparent from the comparisons in Figs. 7-9 that, the radius enrichment strategy yields the most accurate results for a given mesh resolution, with a cost of introducing more extra DOFs (or unknowns in the final linear system) than the two other enrichment strategies. Since all points within a radius R enri to the singular points will be enriched, too many extra DOFs may be introduced if an unnecessarily large R enri has been chosen. For a given R enri , the number of extra DOFs is also larger for a finer mesh. On the other hand, the point enrichment method introduces fewest extra DOFs, but the accuracy is the lowest among the three enrichment methods. From practical application point of view, it is recommended to apply the radius enrichment method with a small enrichment area.
It is of more interest to compare the computational efforts to achieve a similar accuracy. In this regard, we have also plotted in Fig. 11 the L 2 errors of the velocity potential as function of the total number of unknowns, which is an indicator of CPU time. The number of unknowns for different enrichment strategies and different mesh sizes are listed in Table 1 for linear FEMs and Table 2 for quadratic FEMs. It is apparent that the local enrichment increases only marginally the total number of unknowns, while reducing the global errors significantly.
To illustrate how the XFEMs has increased the accuracy of the local flows, the horizontal velocity along the flat plate is shown in Fig. 12 On the contrary, the conventional linear FEM fails to capture the strong variation of the flow at two ends of the plate. Since the applied FEM is only C 0 continuous, the presented velocity in the figure at each point is the average value of the velocities calculated in the elements sharing the point. For the singular elements, the velocity was obtained by differentiating the shape functions in Eq. (27), and we have added more points within the element to better illustrate the variation of velocity therein. In theory, the nodal FEMs based on shape functions are only C 0 continuous, and thus the velocities are discontinuous between elements for both ordinary FEMs and XFEMs. In Fig. 12, the velocity distribution appears smoother for FEM because we have used the average of the velocities evaluated on the adjacent elements as the nodal values. If the averaging is not applied, the velocities will appear discontinuous at all nodes. Since the solution representation is more accurate in the enriched element than that in its neighboring ordinary element, a more obvious jump of the velocity has been observed (at x/a = ±0.75 in the present case).

Heaving rectangular cylinder on free surface
In this part, under the framework of linear potentialflow theory in the frequency domain, a floating heaving rectangular cylinder on a free surface is considered. B and D are used to represent beam and draft of the rectangle, respectively. The considered water depth is h = 40D, and the beam-to-draft ratio B/D is taken as 2.0. An illustration of half of the domain is presented in Fig. 13. In theory, a radiation condition should be applied at x → ±∞. In practice, it is impossible to model an fluid domain with infinite extension, and a truncation at a certain horizontal distance L x from the rectangle must be made. The radiation condition is then applied on a control surface S m , which is chosen sufficiently far from the structure. In this study, we choose a horizontal truncation distance as twice the longest wavelength that will be studied, and use the same computational domain for all different cases.
To reduce the computational costs, the symmetric property of the considered problem is utilized, and thus only half of the fluid domain is considered. At the symmetry plane, horizontal velocity is equal to zero, i.e. ∂φ/∂x = 0. An illustration of the computational domain is presented in Fig. 13.

Linear added mass and damping coefficients
Fig. 14 displays the non-dimensional added mass and damping coefficients for different truncation distances L x from the rectangle. A non-dimensional wave number ω 2 B/(2g) = 0.1 has been considered in the calculations, corresponding to the longest wave that will be considered in this section. If the selected L x has negligible results for the longest wave, it is also considered as sufficient for the shorter waves. It is apparent from the results in Fig. 14 that the hydrodynamic coefficients do not change as long as L x /λ ≥ 1.0. Here λ is the wavelength. L x = 2λ will be applied in our later analyses in this section.
Matched multi-block meshes in the fluid domain are utilized as a starting point, with block I and block II fitted to the body surface, block IV below the body surface, block III and block V away from the structure. See an example of the meshes in Fig. 15, generated from the open source mesh generator GMSH. The following parameters are defined to denote the number of elements along the sides of the blocks to control the mesh densities in different blocks: N rx is the number of elements on the bottom of the rectangle, N ry along the side wall of the rectangle, N ix along the free surface in the inner block, N iy in vertical direction of internal block at symmetry face. Correspondingly, N oy represent the element number in the horizontal direction of external domain on free surface, N oy represent the element number in vertical direction of external domain at symmetry face. Here N rx must be equal to N ix so that blocks I and II will match at their common boundary. For simplicity, we will also take N rx = N ry = N ix = N iy . Meshes in blocks IV and V are stretched along the vertical direction toward the sea bottom using a stretching radio of 1.1. The added mass and damping coefficients are calculated by the four different FEMs, and the results are compared with the experimental results reported in Vugts (1968), the linear numerical potential-flow calculations by Liang et al. (2015). Liang et al. (2015) used the 2D HPC method in the frequency domain, and have taken account of the local singularity by a domain decomposition strategy, where the  Figure 14: Convergence performance of the horizontal length from the rectangle to the matching boundary when the square of forcing frequency is ω 2 B/(2g) = 0.1. A 33 =heave added mass, B 33 =heave radiation damping, S=submerged cross-sectional area, ρ=mass density of water, ω=circular frequency, Lx = horizontal length from the rectangle to the matching boundary, λ = wavelength of radiated waves.
local corner-flow solutions were matched with the outer domain represented by the harmonic-polynomial cells. The mesh parameters used in our FEMs are listed in Table 3, in which N p denotes the number of DOFs (including additional DOFs for XFEM) in the computational domain. The present numerical results agree excellently well with those by Liang et al. (2015). All numerical results seem to deviate from with the experimental results at low frequencies. As it is commented in Vugts (1968), the uncertainties in the experimental results for ω 2 B/(2g) < 0.25 may have been high. For ω 2 B/(2g) ≥ 0.25, the numerical results agree better with the experiments. The small differences may have been contributed by the viscous flow separation at the sharp and other nonlinearities which will occur in reality. From the results in Fig. 16, we may conclude that all the numerical methods in the comparison are be able to accurately predict the linear hydrodynamic coefficients with an affordable effort. It is also observed that the XFEMs do not show clear advantages in the linear hydrodynamic analysis, which is expected as only integrals of velocity potential (multiplied by the normal vector) over the mean wetted body surface are involved in the pressure integration. As seen in the corner flow solution in Sect. 2.3 , the velocity potential is not singular at the corner. However, the fluid velocity close to the sharp corners is singular, which poses great challenges in nonlinear wave loads analysis as will be explained further.

The 2nd order mean vertical force
The calculation of 2nd order wave loads based on pressure integration involves the integration of the quadratic terms of fluid velocities on body surface, which are singular but integrable near the sharp corners. Zhao and Faltinsen (1989) showed that the near-field approach based on direct pressure integration without special consideration of the singularity is very difficult to achieve convergent results, and the approach based on momentum and energy relationship or similar were much more efficient and ro- (2) y /(ρω 2 η 2 3a B), F (2) y = 2nd order mean vertical force, ρ= mass density of water, η 3a = heave amplitude, B= beam, n= enrichment function number. Linear XFEM employs mesh 2 in Table 4, quadratic XFEM employs mesh 2 in Table 5. bust. The later approach often involves integration on a control surface and a free surface confined by control surface and structure surface. Similar conclusions have been obtained later by many others (e.g., see Chen, 2007;Shao, 2019;Cong et al., 2020).
The time averaged 2nd order vertical hydrodynamic force acting on the heaving rectangle by direct pressure integration over the mean wet body surface can be expressed as: where S B0 denotes the wetted mean body surface. η 3 is the heave motion define as η 3 = Re[η 3a e iωt ], with η 3a as the heaving amplitude. n 3 is the vertical component of the normal vector. T is the oscillatory period expressed as T = 2π/ω. φ (1) and φ (2) represent the first and second order velocity potential, respectively. A waterline integral due the fluctuation of waves near the mean water level is neglected as it does not contribute to the vertical loads in this particular case. The time derivatives of first and second order velocity potential equal to zero after time average over one period, and thus Eq.(46) can be simplified as: From a theoretical perspective, the near-field approach, far-field approach and the approaches based on control surfaces should be mathematically equivalent. However, since it is very difficult for the conventional numerical (2) y = F (2) y /(ρω 2 η 2 3a B), F (2) y = 2nd order mean vertical force, ρ= mass density of water, η 3a = heave amplitude, B= beam, r= enrichment radius. Linear XFEM employs mesh 2 in Table 4, quadratic XFEM employs mesh 2 in Table 5.
methods, e.g. FEM, FDM and BEM, to accurately describe the exact fluid velocities close to sharp corners, slow grid-convergences are expected for the near-field approach when it is applied to calculate the 2nd order wave loads. Despite difficult, it is still believed by the authors of this paper that, the strong variation of the local velocities can be captured accurately if an appropriate numerical method is adopted, thus the near-field approach can still be a good option for 2nd wave-load analysis. A good example of such a numerical method is the domain decomposition strategy developed by Liang et al. (2015), where the solutions in the local domain surrounding the sharp corners are represented by the analytical corner-flow solutions. The strategy leads to very accurate and efficient near-field result, but it is not easy to implement for general purposes. The XFEM is a more powerful and general-purpose framework, which allows us to easily and explicitly include, for instance the singular corner-flow solutions as enrichment functions, in the local finite-element approximations. It also inherits other good features of the conventional FEMs, e.g. the use of unstructured meshes.
For the considered rectangle, the interior angle at each corner is γ = 90 • , where γ is the interior angle as illustrated in Fig. 2. Eq. (13) presents all possible fundamental solutions to the corner flows, among which we choose only the first a few as our enrichment function. The first term with j = 1 is ϕ = A 1 r 2 3 cos( 2 3 θ), and the resulting radial velocity ∂ϕ ∂r and circumferential velocity 1 r ∂ϕ ∂θ are in the form of r − 1 3 −singularity as r → 0, which are difficult for any regular functions to achieve good approximation.
In Fig. 17, we compare the non-dimensionalF (2) y when different number of terms from Eq. (13) are included as enrichment functions. As shown in Fig. 17(a), for linear XFEM, convergent result has been achieved for enrichment function number n ≥ 3. For quadratic XFEM, the convergence will be achieved with n ≥ 1, as demonstrated in Fig. 17(b). The reason that a linear XFEM needs more enrichment functions than a quadratic XFEM is as follows: the fundamental solution of a corner flow contains a singular term with j = 1 in Eq. (13) and other higher-order nonsingular terms with j ≥ 2. Those non-singular terms are more accurately captured by the regular quadratic shape functions, and thus it seems to be sufficient for a quadratic XFEM to use only the singular enrichment function from Eq. (13). Based on the discussion above and the numerical observation, only three enrichment functions will be considered in later analyses. Adding unnecessarily too many higher-order terms with j > 3 will pose extra difficulties in numerical integration. On the other hand, the extra DOFs due to enrichment will increase rapidly with the number of enrichment function at each enrichment point. Fig. 18 displays the non-dimensional 2nd order mean vertical forceF (2) y for ω 2 B/(2g) = 1.0 as function of 2R enri /B. The numerical results indicate that, for both linear and quadratic XFEMs, the convergence is achieved when 2R enri /B ≥ 0.2. The results also suggest that it is unnecessary to use a too large enrichment radius, because the results do not seem to improve further as long as R enri is greater than a threshold value of approximately 0.2. On the other hand, larger R enri also means more extra DOFs and unknowns.
In Fig. 19(a), the numerical results ofF (2) y by the linear FEM and the linear XFEM are compared with a reference solution in Liang et al. (2015) based on conservation of fluid momentum (CFM). Direct pressure integration has been applied in the present FEM analyses. Mesh 1, mesh 2 and mesh 3 in the parentheses indicate coarse, medium and fine meshes, respectively. Details of the mesh parameters are shown in Table 4. Apparently, the linear XFEM is more accurate than linear FEM as seen from their comparisons with the CFM results (Liang et al., 2015). Convergent result can be achieved rapidly after refine mesh 1 to mesh 2 via linear XFEM. The unknown numbers (or total DOFs) of mesh 1 and mesh 2 are 28866 and 81421 respectively when linear XFEM is applied. On the con- (2) y = F (2) y /(ρω 2 η 2 3a B), F (2) y = 2nd order mean vertical force, ρ= mass density of water, η 3a = heave amplitude, B= beam, ω= circular frequency. The CFM= conservation of fluid momentum.
trary, the linear FEM has not reached the convergence even with the finest mesh, i.e. mesh 4 with total DOFs of N p = 556146 in Table 4. Note that the total number of unknowns, or total DOFs, are different for a FEM and a XFEM, even though the same mesh is used. This is due to the extra DOFs introduced in the XFEM as a result of local enrichment. For quadratic methods, we also consider three different meshes, i.e. coarse, medium and fine meshes, represented by mesh 1, 2 and 3 in Table 5, respectively. As illustrated by the comparisons in Fig.19(b), the conventional quadratic FEM is not able to reach a convergence even with the finest mesh (mesh 3) with a total DOFs of N p = 406631. On the contrary, the quadratic XFEM results are convergent with the medium mesh (mesh 2, N p = 15416). In fact, results based on the coarse mesh (mesh 1, N p = 9293) are already very close to the reference results. In this coarse mesh resolution, only 4 elements are distributed on half of the rectangle bottom.
Comparing the two XFEMs, quadratic XFEM has shown much faster mesh-convergence rate than linear XFEM. More specifically, convergent results can be reached by quadratic XFEM with less than N p = 15416 DOFs, while its takes N p = 81421 for the linear XFEM. Therefore, the quadratic XFEM is considered as more competitive. From the standpoint of solution enrichment, the quadratic XFEM can be seen as a combination of global and local enrichment, with a global enrichment achieved via higher Lagrange polynomials in regular shape functions, and a local enrichment realized by adding prior knowledge to the local approximation space. The linear XFEM, however, only enriches the solution locally. Therefore, it is generally expected that the quadratic XFEM over-performs the linear XFEM. (2) y = F (2) y /(ρω 2 η 2 3a B), F (2) y = 2nd order mean vertical force, ρ= mass density of water, η 3a = heave amplitude, B= beam, ω= circular frequency. The CFM= conservation of fluid momentum.

Application of unstructured meshes
In the previous subsections, a multi-block structured mesh was adopted for demonstration purpose, and the numerical results based on XFEMs were very encouraging. However, it is well-known that, one of the most powerful property of FEM is that it allows for the use of unstructured mesh without having to modify the numerical code. It is much easier for the unstructured meshes to deal with problems involving complex boundaries. In this subsection, the unstructured mesh will be adopted to investigate the same problem that have been studied in the previous subsection.
An example of the unstructured mesh close the 2D rectangle, generated from the open-source mesh generator GMSH, is shown in Fig. 20. The following parameters are defined to control the number of elements on the fluid boundaries: N rx is the number of elements on the bottom of the rectangle, N ry along the side wall of the rectangle, N f x along the free surface, N sy along the symmetry face, N bx along the bottom of the computational domain and N my along the matching boundary. Furthermore, for both linear and quadratic mesh, the mesh is stretched by a fixed stretching radio of 1.1 along the body boundary, so that the meshes are finer close to the corners. The meshes are also stretched vertically towards the bottom of the fluid and horizontally towards the matching boundary, using stretching factors of 1.08 and 1.05 respectively. The meshes are so adapted that the mesh density is higher close to the body and the free surface.
The 2nd order mean vertical force on the heaving rectangle at free surface is studied again in the frequency domain by using the unstructured mesh and the four FEMs, and the corresponding results for linear FEMs and quadratic FEMs are shown in Fig. 21(a) and Fig. 21(b) respectively. The main parameters of the applied unstructured meshes are summarized in Table 6.
Due to the use of unstructured meshes and stretched grid on the fluid boundaries, it is expected that the required total number of unknowns are much smaller than that of the multi-block structured meshes. This has also been confirmed by our numerical results in Fig. 21(a) and Fig. 21(b). As seen in the figures, to achieve convergent results forF (2) y , it is sufficient to use mesh 1 (total DOFs N p = 40854) and mesh 2 (N p = 5870) in Table 6 for linear XFEM and quadratic XFEM, respectively. On the other hand, as expected, the results of the conventional FEMs are not convergent when the same meshes as the corresponding XFEMs are used. There is one point that must be clarified for the results of the conventional linear FEM and quadratic FEM. In Fig. 21, the linear FEM results appear to be closer to reference results than that of quadratic FEM. This is due to the fact that, mesh 1 as used by linear FEM is much finer than mesh 2 used by quadratic FEM.
In Liang et al. (2015), a modified HPC method based on domain decomposition method was developed to solve the same hydrodynamic problem of the heaving rectangle at free surface in the frequency domain. Corner-flow solutions were used in the inner domain surrounding the  sharp corner, while the outer domain solutions were represented by overlapping harmonic-polynomial cells. The inner and outer domain solutions are matched at their common boundaries. This method was shown to be capable of providing convergent 2nd order mean wave loads by using 80 elements along half of the bottom and in total approximately 352000 unknowns, while our linear and quadratic XFEM models need much fewer unknowns (only 5870 for quadratic XFEM and 40854 for linear XFEM) to achieve equally good results. In a nutshell, the superiority of the present study over Liang et al. (2015) is twofold. Firstly, from implementation point of view, the enrichment strategy based on Partition of Unity in XFEMs to include singular functions near the corner is easier and more flexible. Secondly, the unstructured meshes are allowed in XFEMs, which enables XFEMs to deal with more complex structures, whereas it is expected to be more difficult for the HPC method.

Conclusions
The XFEM is applied as an accurate and efficient tool to solve 2D potential-flow hydrodynamic problems for structures with sharp edges. To demonstrate the advantages of XFEM, four FEM codes, in 2D, including the conventional linear and quadratic FEMs, and the two corresponding XFEMs, are implemented and compared. All of our results have confirmed that the XFEM is a promising framework to deal with potential-flow hydrodynamic problems involving structures with sharp edges. Three different enrichment strategies, including: the point enrichment, patch enrichment and radius enrichment, are also investigated in the study of the uniform flow around an infinite-thin flat plate. The first two enrichment methods are found to be mesh-dependent, and are not able to achieve the expected spatial convergence rate. However, the radius enrichment method is mesh-independent, and has shown remarkably better accuracy and spatial convergence rate. Therefore, it is considered as the best option over the other three counterparts. By studying the horizontal fluid velocity along the flat plate, we also demonstrate that XFEMs are capable of capturing the strong flow variation close to the endpoints, which cannot be represented by the conventional FEMs.
For a heaving rectangular cylinder on the free surface, both the conventional FEMs and XFEMs can accurately predict the linear hydrodynamic coefficients with acceptable computational efforts, indicating that the singularity at sharp corner is inconsequential to the linear hydrodynamic loads. However, it has important effects on the 2nd order mean wave loads if the direct pressure integration is employed, because the singular flow velocities are involved. Compared with reference results based on conservation of fluid momentum, both linear and quadratic XFEMs have shown encouraging results even with a relative coarse mesh resolution, while the quadratic XFEM has an overall better performance than the linear XFEM. On the contrary, it is difficult for the two conventional FEMs to achieve convergence even with an extremely fine mesh. For the quadratic XFEM, it is also found sufficient to only include the first singular term from the corner-flow solutions in the local enrichment, while it is beneficial to include a few more, e.g. 3 terms, in the local enrichment in the linear XFEM.
As a final demonstration, we show that the adoption of unstructured meshes and local refinement close to the sharp edges have a great potential to further reduce the total number of unknowns to achieve a desired accuracy.

Appendix A. shape function and isoparametric element
Commonly, the shape function is defined in an element, for simplicity, using n e p to represent the number of the nodes in a single element. Referring to Zienkiewicz et al. (2005), for 4-node quadrilateral linear FEM, namely n e p = 4, the shape function defined on a parametric ξη-plane can be written as: (1 + ξ i ξ) (1 + η i η) , i = 1, · · · , 4.
where (ξ i , η i ) denote the normalized coordinates at node i.