Upwind finite element-PML approximation of a novel linear potential model for free surface flows produced by a floating rigid body

Applied Mathematical Modelling 103 (2022) 824–853 Contents lists available at ScienceDirect Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm Upwind finite element-PML approximation of a novel linear potential model for free surface flows produced by a floating rigid body A. Bermúdez a , b , d , O. Crego a , ∗, A. Prieto c a Departamento de Matemática Aplicada, Universidade de Santiago de Compostela, Santiago de Compostela E-15782, Spain b Instituto de Matemáticas (IMAT), Universidade de Santiago de Compostela, Santiago de Compostela E-15782, Spain c CITIC, Department of Mathematics, Universidade da Coruña, Spain d Instituto Tecnológico de Matemática Industrial (ITMATI), Santiago de Compostela E-15782, Spain a r t i c l e i n f o Article history: Received 26 March 2021 Revised 2 November 2021 Accepted 8 November 2021 Available online 17 November 2021 Keywords: Linear potential model Free surface Kelvin wakes Perfectly matched layer (PML) Galerkin finite element method Upwinding a b s t r a c t A novel linear potential model is presented to compute free surface flows of incompressible fluids produced by the motion of a floating rigid body in the presence of an underlying non-uniform flow. In particular, the proposed model enables the accurate numerical simulation of the Kelvin wake pattern in a computational domain of reduced size. The governing equations are obtained by using an Arbitrary Lagrangian Eulerian (ALE) formulation which involves the underlying velocity of the fluid around the floating body assuming flat free surface, and a non-dimensional analysis to derive the novel linear system of equations for free surface flows. The discretization of the proposed model is made by a standard Galerkin finite element method, where a SUPG-inspired upwinding strategy has been used in combination with a Perfectly Matched Layer technique which allows truncating the original unbounded fluid domain without introducing spurious reflections in the Kelvin wake pattern. The numerical simulations computed with the proposed approach are compared with the results obtained by the classical linear potential model with uniform underlying flow and also with those from the full incompressible Navier–Stokes equations equipped with the k −ωSST turbulent model. This numerical comparison is discussed in terms of a classical hydrodynamic floating body benchmark involving the Wigley hull. ©2021 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) 1. Introduction The applications of hydrodynamics to naval architecture acquired increasing importance since the 1970s (see for instance, the classical textbook by Newman [1] ). In this context, mathematical models based on the theory of potential flow have been widely utilized for the study of ship motions and the wave patterns generated by floating bodies. The purpose of this work is the derivation and numerical simulation of a novel linear potential model, which governs the propagation of the free surface waves generated by a moving floating body in the presence of an additional underlying fluid stream. ∗Corresponding author. E-mail addresses: [email protected] (A. Bermúdez), oscar[email protected] (O. Crego), andr[email protected] (A. Prieto). https://doi.org/10.1016/j.apm.2021.11.004 0307-904X/© 2021 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) A. Bermúdez, O. Crego and A. Prieto Applied Mathematical Modelling 103 (2022) 824–853 Nomenclature a subscript denotes functions in AL coordinates m subscript denotes functions in material coordinates ˆ denotes functions defined in the PML domain ˜ denotes non-dimensional expressions p = (p 1 , p 2 , p 3 ) material coordinates x = (x 1 , x 2 , x 3 ) Eulerian coordinates z = (z 1 , z 2 , z 3 ) AL coordinates { g 1 , g 2 , g 3 } basis of the reference frame in AL coordinates b external body forces d underlying velocity ggravitational acceleration h d mesh tetrahedron diameter k upwinding function n outward unit normal vector to ∂A t n a outward unit normal vector to ∂ A and ∂ A F n ∞ outward unit normal vector to ∂A ∞ oorigin of the reference frame in Eulerian coordinates ttime u displacement of motion Y f v material velocity in Eulerian coordinates v ap apparent velocity in AL coordinates v r ship velocity in Eulerian coordinates v s stream velocity in Eulerian coordinates B w , L w , D w main dimensions of the hull at midship (beam, waterline length and draft, respectively) Fr Froude number F Y f gradient of Y f in AL coordinates H, L , V typical values of height, length and velocity used to get dimensionless variables I identity tensor P X f , P Y f , P Z f reference mappings of motions X f , Y f and Z f R , R ∞ inner and outer radius of the annular PML domain S , C PML tensors (see Sect. 6.2) T trajectory of motion X f T Cauchy stress tensor X f motion of the fluid Y f motion perturbation defined in AL configuration ηa vertical component of Y f Z f underlying motion βauxiliary variable for upwinding γ, ˆ γ, ˆ r , σθPML functions δSUPG SUPG constant ηfree surface elevation in Eulerian coordinates πpressure in Eulerian coordinates ρf mass density of the fluid σr radial PML absorption profile ϕvelocity potential in Eulerian coordinates φa velocity potential perturbation in AL coordinates (r, θ, z) cylindrical coordinates in AL configuration Domain nomenclature material (or Lagrangian) domain A t fluid domain in Eulerian coordinates f t , r t , b t free surface, wet surface and seabed boundaries in Eulerian coordinates A fluid domain in AL coordinates f a , r a , b a free surface, wet surface and seabed boundaries in AL coordinates A F truncated fluid domain ∂ A f , ∂ A r , ∂ A b free surface, wet surface and seabed boundaries in the truncated domain ∂A i inner boundary between the fluid and the PML domains 825 A. Bermúdez, O. Crego and A. Prieto Applied Mathematical Modelling 103 (2022) 824–853 A ∞ PML domain ∂ A f ∞ , ∂ A b ∞ , ∂ A ∞ free surface, seabed and outer boundaries of the PML domain Recently, the increase of computational resources has allowed to solve numerically more complex hydrodynamic models such as the Reynolds-averaged Navier–Stokes equations (see Kinaci et al. [2] ). However, accurate numerical approximations of these models are computationally demanding. Consequently, simplified models, Kuznetsov et al. [3] , are often used in the early stages of ship design, since they are suitable for quick (although may be less accurate) computations obtained in a feasible time. In this work as in many others, the Navier–Stokes equations are simplified into the Bernoulli and Laplace equations assuming irrotational flow of an incompressible and inviscid fluid. This manuscript presents a rigorous approach using an Arbitrary Lagrangian–Eulerian (ALE) formulation together with a linearization to obtain a novel linear potential model. In the case of unsteady flow phenomena, the description of their governing equations using a Lagrangian or an ALE framework enables to work in a fixed (i.e., time independent) computational domain (see, for instance, some applications to the Navier–Stokes equations in Sarrate et al. [4] , Benítez and Bermúdez [5] ). In previous works, linear models similar to but different from the one proposed in this document (see for instance [6–8] ) have been deduced using formal arguments rather than a rigorous mathematical justification (see Sections 3 and 4 ). An important feature of the faced problem is that the free surface perturbation associated with the body motion can be propagated far from their source (see Fig. 1 ), which can represent a drawback when the original unbounded fluid domain is truncated for discretization purposes. Currently, the Perfectly Matched Layer (PML) is a standard technique to model wave propagation phenomena in unbounded domains, in such a way that the computational domain of interest is surrounded by an absorbing artificial layer, which does not introduce any spurious reflections in the solution of the original problem. It was first proposed by Berenger [9] for electromagnetic problems but it is also widely used in other fields as acoustics (see Abarbanel et al. [10] , Qi and Geers [11] , Bermúdez et al. [12] ), or linearised water waves (see Cohen and Imperiale [13] ). In the present approach the two-dimensional Cartesian PML model has been extended to three dimensions in cylindrical coordinates with the aim of absorbing the outgoing wave pattern generating by the floating body (i.e., the Kelvin wake) and preserving the structure of the convected terms presented in the free boundary condition. Moreover, in the 1980s, the Streamline Upwind Petrov–Galerkin (SUPG) method was introduced in Brooks and Hughes [14] for upwinding the convective term in convection dominated flows. The SUPG method can be understood as a stabilization technique (see, for instance, Giuliani et al. [8] , Pironneau [15] , Mola et al. [16] ) for standard finite element discretizations. In order to deal with the convective term involved in the free surface boundary condition, the SUPG method must be applied both in the fluid and in the PML governing equations without introducing numerical artifacts in the discrete solutions. With this purpose, a novel SUPG-inspired upwind strategy has been considered, ensuring the compatibility and the accuracy of the upwinding approach with the free surface boundary in the fluid and PML domains. The structure of the present manuscript is described as follows: First, the motion equation and the constitutive laws for an incompressible and inviscid fluid are described in Eulerian coordinates in Section 2 . Then, an Arbitrary Lagrangian– Eulerian formulation of the fluid governing equations is written in Section 3 , which is then linearised in Section 4 . Afterwards, the corresponding steady state model is introduced in Section 5 to complete the modelling part of this manuscript. Section 6 is focused on the discrete procedure used to solve the proposed model: a standard piecewise linear finite element method in combination with the above mentioned SUPG-inspired upwinding method and the PML technique. Then, in Section 7 the numerical results obtained with the proposed approach are validated under two different perspectives: First, they are compared with respect to the experimental data available in the classical Wigley hull benchmark (see Kajitani Fig. 1. Kelvin wake pattern generated by a cruise ship. 826 A. Bermúdez, O. Crego and A. Prieto Applied Mathematical Modelling 103 (2022) 824–853 et al. [17] ). Second, the numerical results of the present approach are analysed in comparison with those obtained with a full discretization of the unsteady Navier–Stokes equations equipped with the k −ωSST turbulent model. Finally, Section 8 collects the most relevant conclusions. 2. Model in Eulerian coordinates Firstly, the governing equations of the motion of an incompressible and inviscid fluid with a free surface are described in Eulerian coordinates. With this purpose, some preliminary concepts on continuum mechanics will be recalled. Let X f : ×R + → A t be the motion of a fluid past a floating body such that, for each time t , X f (·, t ) is a deformation of the reference (material) configuration, , into the region of the affine space occupied by the fluid at time t, A t , called spatial or Eulerian configuration at time t(see, for instance, Gurtin [18] 1 ). The trajectory of the motion is the set T := { (x, t) : x ∈ A t , t > 0 } and the so-called reference map of the motion, P X f : T → , is defined by P X f (·, t) := (X(·, t)) −1 . Finally, the Eulerian description of the material velocity, this is, the time derivative of the motion written in Eulerian coordinates is given by v (x, t) := ∂X f ∂t (p, t) p= P X f (x,t) . Under the classical Cauchy’s hypothesis in continuum mechanics, the non-conservative form of the motion equation written in Eulerian coordinates is given by ρf ˙ v = div T + b in A t , t > 0 , where T is the Cauchy stress tensor, b is the external body force density and ρf is the mass density of the fluid, which is assumed constant. Since the fluid is considered incompressible and inviscid, it holds div v = 0 , T = −πI , where πdenotes the pressure field. Besides, the body force will be only driven by the gravity effects, so b (x, t) = grad (−ρf gx 3 ) . Thus, collecting the equations written above, the motion equation becomes ρf ˙ v + grad π= grad (−ρf gx 3 ) in A t , t > 0 . (1) Since all domains A t are assumed to be simply connected, which is a plausible assumption for both floating or submerged body in three spatial dimensions, then the flow is shown to be potential, i.e., there exists a scalar field ϕ : T → R such that v = grad ϕ. Consequently, the incompressibility condition and (1) can be written in terms of the velocity potential: ϕ = 0 in A t , t > 0 , ρf ( grad ϕ ) ·+ grad π= grad (−ρf gx 3 ) in A t , t > 0 . (2) Now, boundary conditions must be introduced to complete the mathematical model. The boundary of the unbounded fluid domain is split in three disjoint parts: the free boundary (on the top) of the fluid domain, f t , the wet boundary of a floating (or submerged) body, r t , and the bottom boundary of the fluid domain, b t . Since the fluid is supposed to be inviscid, the normal velocity has to be continuous with respect to the imposed velocity of the floating body as it is assumed to be impervious. Thus, the boundary condition can be written as ∂ϕ ∂n = v r ·n on r t , t > 0 , (3) where v r denotes the possibly time-dependent velocity field of the rigid body. Analogously, on the rigid bottom surface, b t , it holds ∂ϕ ∂n = 0 on b t , t > 0 . (4) These boundary conditions need to be complemented with an appropriate radiation condition at | x −o| → ∞ . Assuming the velocity field is constant at infinity, the radiation condition imposes that grad ϕ = v s for | x −o| → ∞ , t > 0 , (5) being v s the constant stream velocity and othe origin in the Eulerian coordinate system. The hydrodynamic model written above has to be completed by including two boundary conditions associated with the free boundary, as well as adequate initial conditions. The kinetic condition on the fluid free surface models the balance of forces by means of the Bernoulli equation, whereas the kinematic condition reflects the fact that a particle on the free surface, f t , remains there at any time t > 0 . The precise description of these two boundary conditions is postponed and they will be introduced once the model (2) –(5) has been written in the AL configuration. 1 In [18] an introduction to the continuum mechanics framework can be found; the notations in this book will be followed throughout this article 827 A. Bermúdez, O. Crego and A. Prieto Applied Mathematical Modelling 103 (2022) 824–853 3. Arbitrary Lagrangian–Eulerian model To achieve a more tractable model where the physical domain of the hydrodynamic problem does not depend on time, an ALE formulation is used. With this purpose, the motion of the fluid X f : ×R + → A t is given as the composition of two mappings (see Fig. 2 ), Y f : A ×R + → A t and Z f : ×R + → A , as follows: X f (p, t) = Y f (Z f (p, t) , t) . (6) Thus, on the one hand Z f is identified as an underlying motion that determines the Arbitrary Lagrangian (AL) configuration. More precisely, the Eulerian (with respect to Z f ) domain occupied by the fluid in the underlying motion is fixed in time, namely, Z f (, t) = A , and it will be taken as the AL configuration. Therefore, this assumption implies the fixed free surface of the fluid domain in the AL configuration is not being perturbed by the underlying motion. For the sake of completeness, the ALE model is developed in Section 3.1 . On the other hand, motion Y f is identified as a small perturbation from the underlying motion Z f . Consequently, the Kelvin wake pattern is obtained in the ALE formulation as the free surface elevation associated with Y f . Now, the Eulerian potential model (2) –(5) stated in terms of the motion X f is rewritten in the AL configuration. For this purpose, it is necessary to introduce some notations. The velocity of the underlying flow is denoted by d (p, t) := ∂Z f ∂t (p, t) and F Y f (z, t) := grad a Y f (z, t) denotes the deformation gradient of motion Y f , where grad a (respectively, grad a ) denotes the gradient of a vector field (respectively, a scalar field) defined in the AL configuration with respect to variable z. Besides, for a spatial field f : T → R and a material field h : ×R + → R , the subscript a denotes their respective AL descriptions, namely, f a ( z, t ) = f Y f ( z, t ) , t and h a ( z, t ) = h P Z f ( z, t ) , t = h P X f Y f ( z, t ) , t , t , where P Z f (·, t) := Z f (·, t) −1 is the reference mapping of motion Z f . Analogously, the reference mapping associated with Y f is defined as P Y f (·, t) = Y f (·, t) −1 . A straightforward application of an adequate change of variables on the model (2) –(5) leads to the following result: Lemma 3.1. Let ϕbe a solution of model (2) –(5) , which is the potential velocity field associated with motion X f . If Y f , Z f , and X f are motions satisfying (6) , and F Y f , F Z f and F X f are their corresponding gradients then it holds ( ϕ ) a = div a det (F Y f ) F −1 Y f F −t Y f grad a ϕ a = 0 in A ×R + , ∂Y f ∂t = F −t Y f grad a ϕ a −F Y f d a in A ×R + , F −1 Y f F −t Y f grad a ϕ a ·n a = F −1 Y f v r a ·n a on r a ×R + , F −1 Y f F −t Y f grad a ϕ a ·n a = 0 on b a ×R + , F −t Y f grad ϕ a = v s a for | z −o a | → ∞ , t > 0 . In order to write the boundary conditions imposed on the free boundary of the fluid domain, the balance of pressures and loads should be considered. More precisely, a version of the Bernoulli’s theorem in the AL configuration is obtained (see its proof in Appendix A ). Fig. 2. ALE mapping diagram among the different configurations. X f (·, t) is the total motion at time t . Z f (·, t ) and Y f (·, t ) are respectively the underlying motion and its perturbation at time t. Coloured arrows highlight the underlying velocity field in the AL configuration domain, A , and the velocity field in the spatial configuration domain at time t, A t . 828 A. Bermúdez, O. Crego and A. Prieto Applied Mathematical Modelling 103 (2022) 824–853 Theorem 3.1 (Bernoulli equation in AL coordinates) . Under the same hypothesis of Lemma 3.1 , the following equality holds: −ρf  grad a ∂u ∂t t + ( grad a ( grad a u ) t ) d a  F −t Y f grad a ϕ a + ρf grad a grad a ϕ a d a + grad a ρf ∂ϕ a ∂t + πa + ρf gηa = 0 in A ×R + , where ηa (z, t) = (Y f (z, t) −o a ) ·g 3 is the third coordinate of motion Y f and u (z, t) = Y f (z, t) −zis the displacement field of motion Y f . Thus, collecting the results of Lemma 3.1 and Theorem 3.1 , the hydrodynamic model (2) –(5) written in the AL configuration is given by ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ div a det (F Y f ) F −1 Y f F −t Y f grad a ϕ a = 0 in A ×R + , ∂Y f ∂t = F −t Y f grad a ϕ a −F Y f d a in A ×R + , −ρf  grad a ∂u ∂t t + ( grad a ( grad a u ) t ) d a  F −t Y f grad a ϕ a + ρf grad a grad a ϕ a d a + grad a ρf ∂ϕ a ∂t + πa + ρf gηa = 0 in A ×R + , F −1 Y f F −t Y f grad a ϕ a ·n a = F −1 Y f v r a ·n a on r a ×R + , F −1 Y f F −t Y f grad a ϕ a ·n a = 0 on b a ×R + , F −t Y f grad a ϕ a = v s a for | z −o a | → ∞ , t > 0 . (7) 3.1. The underlying motion The hydrodynamic model (7) written in the AL configuration requires the a priori knowledge of the underlying velocity field d a , which is generated due to the ship motion (with a fixed velocity v r ) and a fluid stream velocity (with a fixed velocity v s ). To compute this underlying flow, the (Eulerian with respect to Z f ) domain with flat surface, A , is considered. As it has been mentioned above, this domain is also assumed as the Arbitrary Lagrangian (AL) configuration for the perturbed flow Y f . Under the same hypothesis of incompressibility and non-viscosity of the fluid, the underlying flow is also assumed potential, this is, d a = grad a ϕ 0 . Since the velocity field d a is tangent to the flat free surface, the potential field satisfies (see, for instance, Landau and Lifshitz [19] ) ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ a ϕ 0 = 0 in A , ∂ϕ 0 ∂n a = 0 on f a ∪ b a , ∂ϕ 0 ∂n a = v r a ·n a on r a , grad a ϕ 0 = v s a for | z −o a | → ∞ , (8) where recall that v r a is the ship velocity and v s a is the stream velocity. Notice that ϕ 0 is the so-called double-body potential (see, for instance, Shahjada Tarafder and Suzuki [6] ). The potential field ϕ 0 is decomposed by using the stream velocity as ϕ 0 = v s a ·(z −o a ) + φ0 . Hence, the model can be rewritten in terms of the potential φ0 , the stream velocity v s a , and the apparent velocity, v ap = v s a −v r a as follows: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ a φ0 = 0 in A , ∂φ0 ∂n a = 0 on f a , ∂φ0 ∂n a = −v s a ·n a on b a , ∂φ0 ∂n a = −v ap ·n a on r a , grad a φ0 = 0 for | z −o a | → ∞ . (9) Consequently, the underlying velocity is d a = grad a ϕ 0 = grad a φ0 + v s a . Notice that problem (9) does not involve any convective term and therefore an accurate discretization of this problem does not require any upwinding method. Moreover, in absence of transport phenomena in this model, the decay of the Laplacian solution φ0 avoids the use of large computational domains. Indeed, it has been numerically verified that truncating the domain without any damping methodology does not introduce significant numerical errors, even in the case where this truncation is done at a small finite distance of the floating body. 829 A. Bermúdez, O. Crego and A. Prieto Applied Mathematical Modelling 103 (2022) 824–853 4. Linear approximation Since the motion X f has been assumed as a small perturbation Y f from the underlying motion Z f , a linear model can be used to approximate the velocity potential ϕ a . Associated with Y f in (7) , the linear model is obtained by applying a dimensional analysis to the equations in (7) and keeping only the first order terms with respect to a small parameter. The typical value of the displacement perturbation (which is closely related with the magnitude of the oscillations of the free surface of the fluid) and the length of the floating body are denoted by Hand L , respectively. Besides, the typical value of the apparent velocity of the underlying motion (more precisely, the subtraction of the stream and the rigid floating body velocities) is denoted by V . Then, L/V is the typical time the fluid takes to travel along the floating body. Accordingly, the following non-dimensional variables are introduced: ˜ t = V L t, ˜ z −o a = 1 L (z −o a ) , ˜ v r a = 1 V v r a , ˜ v s a = 1 V v s a , ˜ ηa = 1 L ηa , ˜ d a = 1 V d a , ˜ u = 1 H u , ˜ ϕ a = 1 V L ϕ a , ˜ πa = L ρf V 2 L πa . Hence, once the second equation of model (7) is rewritten in terms of the displacement u (notice ∂Y f ∂t = ∂u ∂t ) the nondimensional version of (7) is given by ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩  div a det ( ˜ F Y f ) ˜ F −1 Y f ˜ F −t Y f  grad a ˜ ϕ a = 0 in ˜ A ×R + , V H L ∂ ˜ u ∂ ˜ t = V ˜ F −t Y f  grad a ˜ ϕ a −V ˜ F Y f ˜ d a in ˜ A ×R + , −V ρf  V H L 2  grad a ∂ ˜ u ∂ ˜ t t + V H L 2 (  grad a (  grad a ˜ u ) t ) ˜ d a  ˜ F −t Y f  grad a ˜ ϕ a + V 2 L ρf  grad a  grad a ˜ ϕ a ˜ d a + 1 L  grad a V 2 ρf ∂ ˜ ϕ a ∂ ˜ t + V 2 ρf ˜ πa + Lρf g ˜ ηa = 0 in ˜ A ×R + , V ˜ F −1 Y f ˜ F −t Y f  grad a ˜ ϕ a ·˜ n a = V ˜ F −1 Y f ˜ v r a ·˜ n a on ˜ r a ×R + , V ˜ F −1 Y f ˜ F −t Y f  grad a ˜ ϕ a ·˜ n a = 0 on ˜ b a ×R + , V ˜ F −t Y f  grad a ˜ ϕ a = V ˜ v s a for | ˜ z −o a | → ∞ , ˜ t > 0 , (10) where  ·denotes non-dimensional differential operators and fields. Since H L in typical naval architecture applications (see [1] ), it is assumed that = H L is a small parameter. Consequently, the tensor field ˜ F Y f = I + H L  grad a ˜ u and the Bernoulli equation (third equation in (10) ) can be rewritten as ˜ F Y f = I + O () and V 2 L ρf  grad a  grad a ˜ ϕ a ˜ d a +  grad a ρf ∂ ˜ ϕ a ∂ ˜ t + ρf ˜ πa + O () +  grad a ρf g ˜ ηa = 0 , respectively. Thus, the following linear approximations of the equations written above in their dimensional form are deduced: F Y f = I in A ×R + , (11) ρf grad a grad a ϕ a d a + grad a ρf ∂ϕ a ∂t + πa + ρf gηa = 0 in A ×R + . (12) Now, taking into account the linear approximations (11) and (12) , a linearised model can be written from (7) . For this purpose, F Y f is approximated by I in (10) unless in the second term on the right-hand side of the second equation in (7) . Notice that even if F Y f is kept there, the model will be linear. However, if any other occurrence of F Y f were kept then the model would be non-linear. In addition, notice that if F Y f were approximated by I in all occurrences then ϕ a = ϕ 0 so the Kelvin wave pattern would not be a solution of the linear model derived from (10) . As a consequence, the following linear approximate model is finally adopted: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ a ϕ a = 0 in A ×R + , ∂Y f ∂t = grad a ϕ a −grad a Y f d a in A ×R + , ρf grad a grad a ϕ a d a + grad a ρf ∂ϕ a ∂t + πa + ρf gηa = 0 in A ×R + , ∂ϕ a ∂n a = v r a ·n a on r a ×R + , ∂ϕ a ∂n a = 0 on b a ×R + , grad a ϕ a = v s a for | z −o a | → ∞ t > 0 , (13) 830 A. Bermúdez, O. Crego and A. Prieto Applied Mathematical Modelling 103 (2022) 824–853 where the velocity potential ϕ a , pressure πa and the motion Y f are the unknown fields to be computed, whereas d a is assumed to be known (recall that it is computed from the solution φ0 of problem (9) ). Notice that model (13) is overconstrained. In order to get a balanced model and reduce the unknowns only to the velocity potential ϕ a , in the next section, some analytic computations are made previous to its numerical solution. 5. Steady state Now, once the AL description of the time-dependent hydrodynamic problem has been introduced and linearized in the previous section, it will be solved at the steady state, as the main goal of this work consists in the numerical simulation of the Kelvin wakes and this phenomenon is steady in the proposed AL framework. Hence, considering that velocities v s a and v r a are time independent, the steady state model is deduced from (13) by suppressing the terms involving time derivatives. Unlike the classical hydrodynamic model (see Remark 5.1 below), the system of equations (13) also involves the vector unknown Y f and the pressure πa , in addition to the scalar unknown ϕ a . In order to consider only one scalar unknown, the steady version of the second equation in (13) is multiplied by the outward unit normal vector to the free surface in the AL configuration (for simplicity, it is assumed to be the third vector of the standard Cartesian basis in the AL coordinate system, g 3 ). Thus, grad a Y f d a ·g 3 = grad a ϕ a ·g 3 in A and since grad a Y f d a ·g 3 = d a ·( grad a Y f ) t g 3 = d a ·grad a (Y f ·g 3 ) , then grad a ϕ a ·g 3 = grad a ηa ·d a in A . (14) Similarly, the steady version of the third equation of model (13) is scalarly multiplied by d a to get ρf grad a grad a ϕ a d a ·d a + grad a (πa + ρf gηa ) ·d a = 0 in A . Taking into account (14) and the assumption of constant mass density of the fluid, it leads to ρf grad a grad a ϕ a d a ·d a + grad a πa ·d a + ρf g grad a ϕ a ·g 3 = 0 in A . (15) On the free surface f a , pressure πa is supposed to be spatially homogeneous (usually equal to the atmospheric pressure). Since the outward unit normal vector to f a is n a = g 3 and d a is tangent to the free surface then d a ·g 3 = 0. Due to all the arguments described above, grad a πa ·d a = 0 on f a . Therefore, restricted to the free boundary, Eq. (15) yields the boundary condition ∂ϕ a ∂n a = −1 g grad a grad a ϕ a d a ·d a on f a . Keeping in mind the radiation condition, the potential ϕ a is decomposed using the same procedure applied to compute the underlying motion (see the decomposition for ϕ 0 described in Section 3 ) as ϕ a = v s a ·(z −o a ) + φa . Thus, grad a φa = grad a ϕ a −v s a and the normal derivative of φa can be written in terms of ϕ a as follows: ∂φa ∂n a = ∂ϕ a ∂n a −v s a ·n a . Since v s a is tangent to the free surface boundary, then v s a ·n a = 0 on f a . Moreover, as it is spatially constant then ∂φa ∂n a = −1 g grad a grad a φa d a ·d a on f a . Additionally, the equation written above is equivalent to ∂φa ∂n a = −1 g grad a ( grad a φa ·d a ) ·d a + 1 g grad a φa ·grad a d a d a on f a . (16) In conclusion, the AL formulation (13) for the steady-state becomes ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ a φa = 0 in A , ∂φa ∂n a = −1 g grad a ( grad a φa ·d a )·d a + 1 g grad a φa ·grad a d a d a on f a , ∂φa ∂n a = −v ap ·n a on r a , ∂φa ∂n a = −v s a ·n a on b a , grad a φa = 0 for | z −o a | → ∞ , (17) where recall that v ap = v s a −v r a is the apparent velocity of the fluid observed from the floating body. Remark 5.1. Notice that the present model is a generalization of the classical one found in the bibliography (see, for instance, Newman [1] ) which has been widely used to compute Kelvin wake patterns (see, for instance, Giuliani et al. [8] ). 831 A. Bermúdez, O. Crego and A. Prieto Applied Mathematical Modelling 103 (2022) 824–853 The difference with the model introduced in (17) lies in the free surface boundary condition. Indeed, since v ap is constant in these references, replacing d a by v ap in the free boundary condition of (17) leads to ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ a φa = 0 in A ∂φa ∂n a = −1 g grad a ( grad a φa ·v ap )·v ap on f a , ∂φa ∂n a = −v ap ·n a on r a , ∂φa ∂n a = −v s a ·n a on b a , grad a φa = 0 for | z −o a | → ∞ , (18) which is an equivalent writing of the classical linear model. Thus, the difference between (17) and (18) lies on the type of assumed underlying motion: d a = v ap is constant in (18) , but the underlying motion d a used in (17) varies spatially since it comes from solving (9) , this is, d a = grad a φ0 + v s a .  From a weak formulation point of view, model (17) could lead to an ill-posed variational problem stated in standard Sobolev spaces, as it involves a second order differential operator applied on the free surface f a . However, since d a is a tangential vector to f a , the trace grad a φa ·d a can be defined for φa ∈ H 1 (A ) and an analogous argument could be applied to grad a ( grad a φa ·d a )·d a as far as grad a φa ·d a ∈ H 1 (A ) . Moreover, the second term on the right-hand side of the free boundary condition makes sense for φa ∈ H 1 (A ) because grad a d a d a is also a tangential vector to f a . Consequently, to enforce that grad a φa ·d a belongs to H 1 (A ) , a new unknown βis introduced in model (17) , which is defined by β= grad a φa ·d a . Thus, problem (17) is written including the scalar field βas follows: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ a φa = 0 in A , β= grad a φa ·d a in A , ∂φa ∂n a = −1 g grad a β·d a + 1 g grad a φa ·grad a d a d a on f a , ∂φa ∂n a = −v ap ·n a on r a , ∂φa ∂n a = −v s a ·n a on b a , grad a φa = 0 for | z −o a | → ∞ . (19) 6. Numerical resolution In this section the numerical resolution of the steady problem (19) is discussed. With this purpose, the transport term grad a φa ·d a should be discretized adequately and the unbounded domain A truncated taking into account the radiation condition. On the one hand, since problem (19) is discretized by using a standard finite element procedure, the transport term must be conveniently treated by applying an upwinding methodology. On the other hand, a Perfectly Matched Layer technique is introduced to handle the radiation condition and to truncate the otherwise unbounded domain A . In order to discretize the transport term properly, the second equation of (19) will be rewritten. The proposed upwinding strategy is described in Section 6.1 . This strategy is inspired by the Streamline Upwind Petrov–Galerkin (SUPG) method. In fact, it has been derived because the straightforward application of the standard SUPG method produces spurious numerical oscillations when it is combined with the PML technique. Once the bounded computational domain is described after introducing the PML layer in the original problem (19) , the weak form associated with the resulting system will be written. 6.1. Upwind methodology for handling the transport term The SUPG method (see Brooks and Hughes [14] ) can be understood as a replacement of continuous weighting functions with discontinuous ones in the discretization of a Petrov–Galerkin formulation. Instead of working at the discrete finite element level, firstly an equivalent upwinded continuous problem is introduced and then a standard Galerkin method is directly applied. With the aim of describing the upwinding strategy, the second equation of problem (19) is stated in an arbitrary domain D, β= grad a φa ·d a in D. Now, assume for a moment that ϕ a is given in H 1 (D) . The standard weak formulation of this equation is  Find β∈ L 2 (D) such that  D βξdV z =  D grad a φa ·d a ξdV z , ∀ ξ∈ L 2 (D) . (20) 832 A. Bermúdez, O. Crego and A. Prieto Applied Mathematical Modelling 103 (2022) 824–853 Fig. 5. Mean relative difference computed from six different numerical simulations plotted with respect to σ0 (Right: range σ0 ∈ [10 −3 , 2 ×10 −1 ] ; left: range σ0 ∈ [10 −3 , 10 2 ] ). where σ0 is a given positive parameter. In this manner, the pointwise values of the singular PML absorption profile and its polynomial approximation are identical at the quadrature nodes used in the finite element discretization. Moreover, the PML coefficients involved in the weak formulation are well defined on the entire discrete computational domain. In order to select an optimal value for to σ0 in the PML absorption profile, a reference solution was computed in a very large domain, which is truncated with homogeneous boundary conditions. For this purpose, a computational fluid domain with an exterior radius of 40 m has been chosen to compute the reference solution. Then, a quantitative comparison has been made to evaluate the relation between the reference result (computed in the large domain without PML), φref , and the solution φof the problem (19) (computed in a smaller domain truncated with a PML layer). More precisely, the L 2 - relative difference between both potentials is computed on the free surface by means of  φref −φ 2  φref  2 =  ∂A f | φref −φ| 2 dA z 1 2  ∂A f | φref | 2 dA z 1 2 (33) in six different numerical simulations with Froude numbers Fr = 0 . 25 , 0 . 267 , 0 . 289 , 0 . 316 , 0 . 354 , 0 . 408 . After a swept in σ0 from 10 −3 to 10 2 , the value σ0 = 0 . 02 has been chosen since it minimizes the arithmetic mean of the L 2 -relative differences computed from the six numerical approximations (see Fig. 5 ). The minimum error for this optimal value of σ0 is approximately 7% whereas the relative error would reach 16 . 8% if the computational domain were truncated without PML (i.e., considering σ0 = 0 on the same mesh). Additionally, Fig. 6 shows a qualitative comparison between the reference potential and the potential obtained with the PML technique ( σ0 = 0 . 02 ) in the smallest computational domain, for the case of Fr = 0 . 25 (notice that only a portion of the largest domain used for computing the reference solution is depicted). The spatial pattern of the absolute difference | φre f − φa | between the reference potential and the approximated potential computed with the PML technique is also reported in the bottom panel of Fig. 6 . 7.3. Numerical comparison with experimental data and full Navier–Stokes approximations In this section, the proposed model, the classical linear potential model (18) , and the full Navier–Stokes with two species are compared with the available experimental data at Fr = 0 . 267 . More precisely, the classical linear potential model has been already described in Remark 5.1 and a turbulence k −ωmodel coupled with the full Navier–Stokes equations (see Wilcox [24] ) was solved in a computational domain with two species: air over the water to model the free surface boundary (using a volume of fluid method, see Hirt and Nichols [25] ). Both fluids have been considered as incompressible. The use of a turbulence k −ωmodel makes it necessary to add boundary conditions for turbulence. For this purpose, the turbulence intensity and the turbulence length scale were set at 2 . 45% and 0.095 m, respectively, on inlet and outlet boundaries. Regarding the used computational domains, the classical linear model and the one proposed in the present work have been discretized using the same mesh. Analogous upwinding and PML techniques have been used in both cases. In the case of the Navier–Stokes model, a prism was chosen, no PML layers are involved and in order to deal with the displacement of the floating body the apparent constant velocity was used to include the forwarding constant velocity of the ship in the stream constant velocity. Analogous to the linear potential models, the symmetry with respect to the longitudinal plane of the Wigley hull has been used to restrict the computational domain to one half. In order to recover the interface between the two fluid phases, a refinement was imposed on the plane z = 0 and, thus, its discretization involves a mesh with 33 millions of tetrahedra approximately. The full two-phase Navier–Stokes model has been solved numerically using Ansys Fluent. A set of control points located on the wake support has been used to define the convergence stop criterion. Once constant values are reached at these points, the solution is assumed to be steady-state and the iterative process is stopped. 839 A. Bermúdez, O. Crego and A. Prieto Applied Mathematical Modelling 103 (2022) 824–853 Fig. 6. Left: potential field, φa , for the experiment of the Wigley hull fixed to sink and trim with Fr = 0 . 25 , obtained using the present model with PML. Right: a zoom of the reference potential field without PML in a big domain with 40 m radius. The maximum number of iterations has been fixed to 10 0 0 0. Regarding the computational cost (CPU time), this numerical simulation needed 3 days of CPU time running on 30 threads in parallel. The computational CPU time has been measured using 4 nodes of a cluster server, where each node has 4 processors Intel(R) Xeon(R) Gold 6126 CPU @ 2.60 GHz and 380 GB of RAM. Fig. 7 shows the numerical approximation of the non-dimensional free surface elevation as a function of the nondimensional position along the Wigley hull waterline, x/L w , for the three models as well as the experimental results, with Fr = 0 . 267 . For the linear models, the finest mesh (shown in Fig. 4 ) has been used and the numerical free surface elevation was computed by ηa = −1 g grad a φa ·v ap , (34) as suggested in Giuliani et al. [8] . In all cases, for a direct comparison with the available experimental data, Kajitani et al. [17] , the non-dimensional quantity ˜ η= 2 gηa / | v ap | 2 has been plotted. In order to compare the flow downstream, Fig. 8 shows the three components of the velocity computed using these three different models. The full two-phase Navier–Stokes model and the proposed model show a similar pattern, except for the region of the wake aligned with the edge on the stern. However, the full two-phase Navier–Stokes model seems to damp the wake oscillations more than the proposed model as it can be observed in Fig. 7 . Regarding the classical linear model (see Fig. 8 ), the wake angle (angle between the main Wigley axis and the line where the maxima/minima of the wake pattern are located) and the amplitude of the velocity potential are smaller than in the other two models. 840 A. Bermúdez, O. Crego and A. Prieto Applied Mathematical Modelling 103 (2022) 824–853 Fig. 7. Non-dimensional free surface elevation associated with the Kelvin wake generated by the Wigley hull with Fr = 0 . 267 (experimental results from Kajitani et al. [17] are depicted with orange dots, classical linear model with a red continuous line, Navier–Stokes results with a green dashed line, and results from the present model with a continuous blue line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 7.4. Numerical comparison with respect to the classical model In Figs. 9–11 , the free surface elevations (computed as it is described in Section 7.3 ) for the six experiments reported in Kajitani et al. [17] are compared with the numerical results of the proposed model and the classical linear model, (18) , for three different meshes using linear finite elements. In all these numerical results, the free surface elevation is computed by means of the gradient of a continuous piecewise linear discretization of the velocity potential φa (see (34) ), what would lead to a discontinuous piecewise constant approximation of ηa . To avoid these discontinuities per element, this gradient discretization is projected into a continuous piecewise linear discrete space, which is used to interpolate the free surface elevation at some points on the Wigley hull waterline. In view of the numerical results reported in Figs. 9–11 , the proposed model exhibits a robust behaviour to compute accurately the free surface elevation along the waterline for the six experiments (independently of the refinement type used in the mesh), while the classical linear model presents instabilities on the stern for some experiments in the coarse mesh. Additionally, the resistance coefficients have been computed using for different meshes the classical, (18) , and the proposed linear model. In all cases the resistance coefficient is given by C PR = R p 1 2 ρf | v ap | 2 S , (35) where Sis the wetted surface area at rest defined by S := 0 . 661 L w (2 D w + B w ) (see Kajitani et al. [17] ). Recall that B w , L w and D w are the main dimensions of the hull (see Section 7.1 ), and the pressure resistance R p is computed by R p =  r t Tn ·−v s | v s | dA x =  r t πn ·v s | v s | dA x =  r a πa det F Y f F −t Y f n a ·e 1 dA z . The ALE pressure field πa is computed from the Bernoulli’s equation written in Eulerian coordinates, as πa = (π) a = ρf 2 | v s | 2 −ρf 2 | grad ϕ| 2 −ρf gx 3 a = ρf 2 | v s a | 2 −ρf 2 | F −t Y f grad a ϕ a | 2 −ρf gηa . Next, F Y f is approximated by the identity and ηa by z 3 leading to the formula, R p ≈ r a ρf 2 | v s a | 2 −ρf 2 | grad a ϕ a | 2 −ρf gz 3 n a ·e 1 dA z . 841 A. Bermúdez, O. Crego and A. Prieto Applied Mathematical Modelling 103 (2022) 824–853 Fig. 8. Numerical approximation of the velocity field for the Kelvin wake generated by the Wigley hull with Fr = 0 . 267 using the full Navier–Stokes model, the present model and the classical linear model (each row uses the same colour scale). Table 1 Resistance coefficient (defined by (35) ) computed with different models, meshes and Froude numbers using linear finite elements. The last row shows the experimental results for comparison purposes. Model Mesh Fr 0.25 0.267 0.289 0.316 Present Coarse mesh 6 . 997 ×10 −4 6 . 979 ×10 −4 1 . 046 ×10 −3 1 . 265 ×10 −3 Present Fine mesh 7 . 178 ×10 −4 6 . 705 ×10 −4 1 . 082 ×10 −3 1 . 258 ×10 −3 Present Very fine mesh 8 . 003 ×10 −4 6 . 722 ×10 −4 1 . 169 ×10 −3 1 . 373 ×10 −3 Classical Coarse mesh 7 . 844 ×10 −4 8 . 525 ×10 −4 1 . 176 ×10 −3 1 . 376 ×10 −3 Classical Fine mesh 7 . 632 ×10 −4 8 . 186 ×10 −4 1 . 228 ×10 −3 1 . 416 ×10 −3 Classical Very fine mesh 9 . 568 ×10 −4 9 . 141 ×10 −4 1 . 360 ×10 −3 1 . 4 4 4 ×10 −3 Experimental results 9 . 41 ×10 −4 8 . 27 ×10 −4 1 . 221 ×10 −3 1 . 789 ×10 −3 This approximation could be the cause of the observed underestimation of the resistance coefficient. Table 1 shows the numerical approximation of the resistance coefficients, which are again compared with the experimental data reported in Kajitani et al. [17] . The resistance coefficient obtained from the full two-phase Navier–Stokes model was 1 . 054 ×10 −3 in the case shown in Section 7.3 ( Fr = 0 . 267 ). Finally, the computational cost and the main characteristics of the three different meshes are shown in Table 2 . The CPU time is averaged among the six experiments described above with different Froude numbers using a computer with 4 processors Intel(R) Xeon(R) Gold 6126 CPU @ 2.60 GHz and 380 GB of RAM. The numerical approximations have been computed sequentially and the linear discrete system has been calculated using a direct solver (MUMPS). 842 A. Bermúdez, O. Crego and A. Prieto Applied Mathematical Modelling 103 (2022) 824–853 Fig. 9. Non-dimensional free surface elevation results at different Froude numbers ( Fr = 0 . 25 , 0 . 267 ) using the meshes shown in Fig. 4 and linear finite elements. They are compared with the experimental data (experimental results from [17] depicted with orange dots, classical linear model with lines in reddish colours and results from the present model with lines in bluish colours). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 843 A. Bermúdez, O. Crego and A. Prieto Applied Mathematical Modelling 103 (2022) 824–853 Fig. 10. Non-dimensional free surface elevation results at different Froude numbers ( Fr = 0 . 289 , 0 . 316 ) using the meshes shown in Fig. 4 and linear finite elements. They are compared with the experimental data (experimental results from Kajitani et al. [17] depicted with orange dots, classical linear model with lines in reddish colours and results from the present model with lines in bluish colours). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 844 A. Bermúdez, O. Crego and A. Prieto Applied Mathematical Modelling 103 (2022) 824–853 Fig. 11. Non-dimensional free surface elevation results at different Froude numbers ( Fr = 0 . 354 , 0 . 408 ) using the meshes shown in Fig. 4 and linear finite elements. They are compared with the experimental data (experimental results from Kajitani et al. [17] depicted with orange dots, classical linear model with lines in reddish colours and results from the present model with lines in bluish colours). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 845 A. Bermúdez, O. Crego and A. Prieto Applied Mathematical Modelling 103 (2022) 824–853 Fig. 12. Non-dimensional free surface elevation results at different Froude numbers ( Fr = 0 . 25 , 0 . 267 ) using the “coarse” mesh and the “fine” mesh shown in Fig. 4 , and quadratic finite elements. They are compared with the experimental data (experimental results from Kajitani et al. [17] depicted with orange dots, classical linear model with lines in reddish colours and results from the present model with lines in bluish colours). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 846 A. Bermúdez, O. Crego and A. Prieto Applied Mathematical Modelling 103 (2022) 824–853 Table 2 Computational cost (CPU time) and mesh characteristics. The second column shows the average CPU time (among the six experiments described above with different Froude numbers). The first three meshes include the PML region but the reference mesh does not. The third and fourth columns show the number of vertices along the length and height of the wet surface of the Wigley, respectively. The fifth column shows the number of nodes (at midship) in the radial wake direction from the stern to the exterior boundary of the computational domain. The last column shows the number of tetrahedra of each mesh. Mesh Average time (s) Vertices Tetrahedra Length Height Wake Coarse 78 148 12 34 0 . 362 ×10 6 Fine 951 296 24 68 2 . 898 ×10 6 Very Fine 752 592 48 184 3 . 710 ×10 6 Reference 3153 296 24 586 1 . 771 ×10 7 Table 3 Resistance coefficient (defined by (35) ) computed with different models, meshes and Froude numbers using quadratic finite elements. The last row shows the experimental results for comparison purposes. Model Mesh Fr 0.25 0.267 0.289 0.316 Present Coarse mesh 7 . 535 ×10 −4 6 . 761 ×10 −4 1 . 171 ×10 −3 1 . 332 ×10 −3 Present Fine mesh 7 . 834 ×10 −4 6 . 852 ×10 −4 1 . 202 ×10 −3 1 . 363 ×10 −3 Classical Coarse mesh 9 . 210 ×10 −4 9 . 376 ×10 −4 1 . 355 ×10 −3 1 . 360 ×10 −3 Classical Fine mesh 9 . 511 ×10 −4 9 . 410 ×10 −4 1 . 409 ×10 −3 1 . 465 ×10 −3 Experimental results 9 . 41 ×10 −4 8 . 27 ×10 −4 1 . 221 ×10 −3 1 . 789 ×10 −3 The last columns in Table 2 contain information about the vertices distribution on the meshes (vertices on the length of the ship, nodes on the height of the ship, and nodes along the radial direction of the wake), and the number of tetrahedra. 7.4.1. Quadratic FE To evaluate the sensitivity of the numerical results with respect to the finite element discretization, the velocity potential has been approximated based on quadratic finite elements. These numerical results have been compared with those obtained with linear finite elements. Figs. 12–14 show the free surface elevations (computed following the procedure described in Section 7.3 ) replicating the six experiments reported in Kajitani et al. [17] . Linear and quadratic approximations are compared using the proposed and the classical linear model, (18) , for three different meshes. From this comparison, it can be concluded the use of a quadratic approximation smooths the high frequency spurious oscillations produced by the projection procedure used to compute the free elevation in the linear finite element approximation (described in Section 7.4 ). Finally, Table 3 shows the resistance coefficient, computed using (35) , for the three different meshes and the classical and the proposed linear model. 8. Conclusions A novel linear potential model for the free surface flow produced by the presence of a floating rigid body in motion have been derived by using an ALE formulation and a subsequent linearisation. The difficulties associated with its numerical resolution have been pointed out and discussed in detail: the convective term stated on the free surface boundary and the truncation of the unbounded physical domain where an outgoing radiation condition is imposed. In order to overcome these difficulties, a new approach to upwind the convective term have been developed, which avoids a spurious numerical behaviour when it is combined with a PML technique. The use of a PML layer with a singular absorbing profile guarantee an accurate damping of the solution in the radial direction even when the computational domain is truncated near the floating body. The proposed methodology has been validated by comparing its numerical results with the available experimental data for the Wigley hull benchmark. Currently, further work is being focused on the extension of the proposed numerical methodology based on a finite element discretization combined with a PML technique to solve numerically the unsteady problem, not only to simulate the Kelvin wake pattern but also to include incoming sea surface waves coupled with the rigid motion of the floating body. 847 A. Bermúdez, O. Crego and A. Prieto Applied Mathematical Modelling 103 (2022) 824–853 Fig. 13. Non-dimensional free surface elevation results at different Froude numbers ( Fr = 0 . 289 , 0 . 316 ) using the “coarse” mesh and the “fine” mesh shown in Fig. 4 , and quadratic finite elements. They are compared with the experimental data (experimental results from Kajitani et al. [17] depicted with orange dots, classical linear model with lines in reddish colours and results from the present model with lines in bluish colours). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 848