- Email: [email protected]
The derivation and the computation of kinematic boundary condition
Mathematics and Computers in Simulation 71 (2006) 62–72
The derivation and the computation of kinematic boundary condition S.H. Leea,b,∗ , B.K. Sonia,b a
Hydrodynamics Research Department, Hyundai Heavy Industries, Co., Ltd., Ulsan, Korea b Mechanical Engineering, University of Alabama at Birmingham, USA
Received 26 July 2005; received in revised form 27 October 2005; accepted 16 November 2005 Available online 27 December 2005
Abstract The generation of waves by an advancing ship is studied using a moving grid approach. The kinematic boundary condition in the physical curvilinear coordinate system is derived by allowing the particle on the upper boundary to move along the generated curved coordinate lines of the flow field. Thus, the grid is moved along the hull surface. An approach based on the idea of the modified partial differential equation was proposed to discretize a hyperbolic partial differential equation. The computational results obtained using a standard finite analytic scheme were compared with the results obtained using a higher order finite analytic scheme. Convergent steady-state wave profiles for the Wigley hull were compared with experimental measurements. © 2005 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: A moving grid approach; Kinematic boundary condition; Hyperbolic partial differential equation; Higher order finite analytic scheme
1. Introduction In incompressible viscous fluid flows, density is constant and pressure is considered to be a primitive variable. The difficulty in simulating incompressible flow fields lies in solving the coupling of velocity and pressure. Free surface prediction around the ship hull and achieving far-field predictions of sufficient accuracy are still a difficult task [12]. In the proceedings of the 23rd international towing tank conference (ITTC), it was reported that free surface fitting techniques can potentially be quite accurate but are less suitable for large free surface distortions, or for situations in which the grid has to be moved along walls of a complicated shape [9]. These illustrate the basic motivation to resolve the problems in numerical ship hydrodynamics. There are two different methods of treating free surfaces, which depend on the grid type: either fixed or moving. In this work, an approach based on a moving grid was applied to solving the free surface flow. The kinematic boundary condition was then solved with the momentum equations in a coupled manner. A three-dimensional incompressible Reynolds-averaged Navier–Stokes (RANS) solver and an unsteady structured grid generation solver were used in this work [13,14]. A new computational code for the kinematic boundary condition is developed. The code is added to the above flow solver. Three-dimensional moving free surface turbulent flow around the Wigley hull is computed.
∗
Corresponding author. Tel.: +82 52 230 5578; fax: +82 52 230 3410. E-mail address: [email protected] (S.H. Lee)
0378-4754/$32.00 © 2005 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2005.11.011
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
63
The kinematic boundary condition was used to update the unknown location of the free surface. The equation of the free surface F (x, y, z, t) was decomposed into a z-coordinate and G(x, y, t). Particles on the upper boundary of the free surface are assumed to move up and down vertically. This assumption is widely used in wave simulation [8,17,11,22,4]. In the computation of free surface flow around an actual ship hull, the geometry of the hull has a fully complex shape near the fore and stern. The stern region near the upper surface has a steep grade. The grid in this region may have a very large skew, and may make it difficult to obtain a reasonable solution with this assumption. To avoid this problem, this work newly derives the kinematic boundary condition in the physical curvilinear coordinate system by allowing the particle on the upper boundary to move along the generated curved coordinate lines in the flow field. In the physical curvilinear coordinate system, each component of the velocity vector has the same direction as that of each curvilinear coordinate line and has a physical value. The physical curvilinear component form of the velocity was first introduced by Truesdell [20]. Lee [13] and Warsi [21] derived formulas needed in order to transform differential equations into a physical curvilinear coordinate system. In this work, Lee [13] and Warsi [21]’s method is used. Another important issue is the discretization of the kinematic boundary condition. This condition is a hyperbolic partial differential equation (PDE). Several approaches have been tried to obtain a reasonable discretization. Tahara and Stern [18] used a backward difference for the x-derivative term and a central difference scheme for the y-derivative term. Later, Tahara and Stern [19] used a Beam and Warming linear-multi-step method. They added implicit and explicit fourth-order artificial dissipations to a difference scheme. Beddhu et al. [1] used an explicit scheme to solve this condition without adding an artificial dissipative term. Farmer et al. [6] used a local time stepping scheme and added background dissipation to prevent decoupling of the solution. Later, Farmer et al. [7] used symmetric-limited positive and upstream-limited positive schemes to avoid overturning of the waves. In this work, an approach based on analytic methods is used to discretize the kinematic boundary condition, which requires no artificial dissipative terms. Recently, Lee and Soni [15] derived higher-order finite analytic schemes based on one-dimensional local analytic solutions to solve three-dimensional PDEs. In the present research, the same higher order schemes were derived and applied to two-dimensional hyperbolic PDEs using the idea of the modified PDE. 2. Physical curvilinear coordinate system In the following analysis, xk are Cartesian coordinates, ξ i are curvilinear coordinates and ξ (i) are physical curvilinear coordinates. The relationship between the coordinates ξ i and ξ (i) is defined as follows: √ ξ (i) = ξ (i) (ξ i ) where ξ (i) = gii ξ i , ξ 1 = ξ 2 = ξ 3 = 1 (1) The coordinates ξ (i) resemble the coordinate stretching in each direction of ξ i . In order to transform the coordinates, the vector dr can be written as: dr = ai dξ i = a(i) dξ (i)
where ai =
∂r ∂r and a(i) = (i) ∂ξ i ∂ξ
(2)
ai are covariant base vectors and a(i) are covariant base vectors in the physical curvilinear coordinate system. The relationship between the two covariant base vectors is obtained from (2) as: 1 a(i) = √ ai gii
(3)
The physical curvilinear components of a velocity vector, u(i) , are defined as the magnitude of the ith component projected onto the ith physical curvilinear coordinate direction: u(i) = u · a(i)
(4)
This form was introduced to give physical meaning to the contravariant component of a velocity lost by full transformation. The relationship between the velocity components in the two different coordinate systems is defined as: u(i) = u(j) a(j) · ei = u(j)
∂xi ∂xi 1 = √ u(j) j ∂ξ (j) gjj ∂ξ
(5)
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
64
u(j) = u(i)ei · a(j) = u(i)
∂ξ (j) ∂ξ j √ = gjj u(i) ∂xi ∂xi
(6)
The velocity components in u(i) in Eq. (5) represent (u, v, w). The repeated indices above imply a sum. Using the general transformation laws, for a scalar φ, the gradient is written as: ∇φ =
∂φ (i) ∂φ a = g(ik) (k) a(i) ∂ξ (i) ∂ξ
where a(i) = g(ik) a(k)
and a(i) =
(7)
√ gii ai
ai are contravariant base vectors and a(i) are contravariant base vectors in the physical curvilinear coordinate system. The relationship between the two contravariant metric tensors in (7) is written as: g(ij) =
√ √ gii gjj gij
where gij =
∂ξ i ∂ξ j ∂xk ∂xk
(8)
Also, the relationship between two covariant metric tensors is written as: g(ij) = √
1 √
gii gjj
gij
where gij =
∂xk ∂xk ∂ξ i ∂ξ j
(9)
g(ij) and g(ij) above are the covariant and contravariant metric tensors in the physical curvilinear coordinate system, respectively. The dot product of the two vectors u and ∇φ is expressed as: u · ∇φ = u(i) a(i) ·
∂φ (j) ∂φ a = u(i) (i) ∂ξ (j) ∂ξ
(10)
3. Kinematic boundary condition The concept of the material surface has been used to derive the kinematic boundary condition, which expresses the fact that the free surface moves with the fluid and deforms in shape as the flow progresses. There is no mass transport through the surface, so its boundary is always composed of the same fluid particles. If the equation of the surface is represented by F = 0, then the equation of the surface at t + t is expressed as: F (r + δr, t + t) = 0
(11)
Using a Taylor’ series expansion, the kinematic boundary condition can be derived by taking the first order components in the vector form: ∂F + u · ∇F = 0 (12) ∂t It is assumed here that flow particles move up and down vertically. If the equation of the surface F (x, y, z, t) is decomposed into the z-coordinate and G(x, y, t), (12) can be written in the Cartesian coordinates system by substituting the z-coordinate and G for the surface, F as follows: ∂G ∂G ∂G +u +v =w ∂t ∂x ∂y
where F (x, y, z, t) = z − G(x, y, t)
(13)
(12) can also be expressed in the physical curvilinear coordinate system. The vector form of (12) is transformed by using (10) as follows: ∂F ∂F + u(i) (i) = 0 ∂t ∂ξ
(14)
If we assume that the flow particles move only vertically, (14) can be rewritten as follows: ∂G ∂G ∂G + u(1) (1) + u(2) (2) = w ∂t ∂ξ ∂ξ
where w = u(i)
∂z ∂ξ (i)
(15)
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
65
Fig. 1. The particles on the upper boundary moving up and down vertically.
Fig. 2. The particles on the upper boundary moving along the generated curved coordinate lines.
Here, the surface equation F was decomposed into its z-coordinate and G(ξ (1) , ξ (2) , t). (15) can be used where the vertical movement is dominant [5]. Fig. 1 shows the particles on the upper boundary moving up and down vertically. When flow particles move along the generated curved coordinate lines ξ (3) , (14) can be rewritten as follows: ∂G ∂G ∂G ∂H ∂H + u(1) (1) + u(2) (2) = u(1) (1) + u(2) (2) + u(3) ∂t ∂ξ ∂ξ ∂ξ ∂ξ
(16)
The surface equation F was decomposed into H(ξ (3) ) and G(ξ (1) , ξ (2) , t). H represents the line in the ξ (3) -coordinate direction. Fig. 2 shows the particles on the upper boundary moving along the generated curved coordinate lines. 4. Discretization methods 4.1. Finite analytic scheme The kinematic boundary condition is coupled with momentum equations. This boundary condition is required to regenerate the grid in the flow fields at every time step. Thus, it is necessary to solve this condition accurately. The scheme used here can be applied to the discretization of hyperbolic PDEs. (16) are rearranged to make the standard form of the finite analytic method developed by Chen and Chen [3]. 0=a
∂G ∂G ∂G + b (2) + c + S ∂ξ (1) ∂ξ ∂t
(17)
where S = −u(1)
∂H ∂H − u(2) (2) − u(3) , ∂ξ (1) ∂ξ
a = u(1) , b = u(2) , c = 1
The second derivative terms are added to the right hand side of (17). To minimize the effects of the second derivative terms, each second derivative term can be multiplied by a very small value. Thus modified, the equation for the kinematic boundary condition can be written as follows: ∂2 G ∂ξ
(1) 2
+
∂2 G ∂ξ
(2) 2
+
∂2 G ∂G ∂G ∂G = A (1) + B (2) + C +S ∂ξ ∂ξ ∂t ∂t 2
(18)
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
66
where
1 (1) ∂H (2) ∂H (3) S= −u , +u −u α ∂ξ (1) ∂ξ (2)
A=
u(1) , α
B=
u(2) , α
C=
1 α
(18) is the standard form of the finite analytic method that solves the elliptic PDE. A and B are assumed to be constant over the local element. A six-point finite analytic scheme based on a one-dimensional analytic solution is derived at the vertex of a cell. (18) is rearranged by decomposing into one-dimensional equations as shown below: ∂2 G 2 ∂ξ (1)
∂2 G 2 ∂ξ (2)
−A
∂G = S1 ∂ξ (1)
(19)
−B
∂G = S2 ∂ξ (2)
(20)
∂G ∂2 G −C = S3 2 ∂t ∂t
(21)
S1 + S2 + S3 = S
(22)
where S1 = − S3 = −
∂2 G ∂ξ (2)
2
∂2 G ∂ξ
(1) 2
+B
∂G ∂2 G ∂G − +C + S, 2 (2) ∂ξ ∂t ∂t
+A
∂G ∂2 G ∂G − + B (2) + S 2 (1) (2) ∂ξ ∂ξ ∂ξ
S2 = −
∂2 G ∂ξ (1)
2
+A
∂G ∂2 G ∂G − +C + S, 2 (1) ∂ξ ∂t ∂t
The source terms S1 , S2 and S3 are assumed to be constant in a local cell. The analytic solution for one-dimensional equation [2] is defined in the present coordinate system using three nodal points. The solution for (19)–(21) is expressed as: (1)
− 1) + b1 −
S1 A
(23)
(2)
− 1) + b2 −
S2 B
(24)
G(ξ (1) ) = a1 (eAξ G(ξ (2) ) = a2 (eBξ
G(t) = a3 (eCt − 1) + b3 −
S3 C
(25)
The coefficients a1 , b1 , a2 , b2 , a3 and b3 are determined by applying the boundary conditions at ξ (1) = hD and −hU , ξ (2) = hN and −hS , respectively. The subscripts U and D denote upstream and downstream, respectively. Fig. 3 shows local elements. The finite analytic solution for the one-dimensional modified PDE can be obtained by taking sufficiently large values of A, B and C as: GP = aD GD + aU GU − CPX S1 If A > 0,
aD = 0, aU = 1, CPX =
(26) hU A
else
aD = 1, aU = 0, CPX =
−hD A
GP = aN GN + aS GS − CPY S2 If B > 0,
aN = 0, aS = 1, CPY =
GP = aT Gn+1 + aB GnP − CPT S3 P aT = 0, aB = 1, CPT =
t C
(27) hS B
else
aN = 1, aS = 0, CPY =
−hN B (28)
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
67
Fig. 3. Local elements.
A six-point finite analytic scheme is derived by substituting the equations for S1 , S2 and S3 from (26)–(28) into (22) as: 1 aD aU aN aS n GP = CP G + GD + GU + GN + GS − S (29) CPT P CPX CPX CPY CPY where 1 1 1 1 = + + CP CPX CPY CPT 4.2. Higher-order finite analytic scheme The one-dimensional equation in (19)–(21) can be solved analytically in a local cell by using one center point and two end cell faces as: Gp = ad Gd + au Gu − Cpx S1
(30)
Gp = an Gn + as Gs − Cpy S2
(31)
Gp = at Gn+1/2 + ab Gn−1/2 − Cpt S3 p p
(32)
ad = at =
1 − e−Ahu , eAhd − e−Ahu 1 − e−Ct/2 , − e−Ct/2
eCt/2
au =
eAhd − 1 , eAhd − e−Ahu
aw =
an =
eCt/2 − 1 , − e−Ct/2
eCt/2
1 − e−Bhs , eBhn − e−Bhs Cpx =
as =
eBhn − 1 , eBhn − e−Bhs
1 hd (e−Ahu − 1) + hu (eAhd − 1) , A eAhd − e−Ahu
1 hn (e−Bhs − 1) + hs (eBhn − 1) t (e−Ct/2 − 1) + (eCt/2 − 1) , C = pt B eBhn − e−Bhs 2C eCt/2 − e−Ct/2 Here, the subscripts u, d, n and s on G represent the values on the cell faces. The value at point p is obtained from the values at the cell faces d, u, n and s. Fig. 3 shows the grid size of a local cell. For sufficiently large values of A, B and C, the analytic coefficients in (30)–(32) can be approximated as: Cpy =
hu −hd else ad = 1, au = 0, Cpx = , 2A 2A hs −hn t = else an = 1, as = 0, Cpy = , at = 0, ab = 1, Cpt = 2B 2B 2C
If A > 0,
ad = 0, au = 1, Cpx =
If B > 0,
an = 0, as = 1, Cpy
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
68
A formula for an unsteady two-dimensional equation is obtained as: ab n−1/2 ad au an as Gp = Cp Gp + Gd + Gu + Gn + Gs − S Cpt Cpx Cpx Cpy Cpy
(33)
Here, the value on cell face u in (33) can be solved analytically for a local cell by applying the boundary condition at points U and D, or UU and P. However, the value G− u on cell face u obtained by applying the boundary condition at points U and D is different from the value G+ u obtained at points UU and P. Therefore, in this work, G on the cell face in (33) will be defined as [15]: Gd =
− G+ d + Gd , 2
Gn−1/2 = p
n−1/2 + )
(Gp
Gu =
− G+ u + Gu , 2
Gn =
− G+ n + Gn , 2
Gs =
− G+ s + Gs 2
n−1/2 − )
+ (Gp 2
(34) (35)
− The analytic solutions for G+ u and Gu in (34) are obtained by applying the boundary condition at GUU and GU , GU and GD , respectively: + + + G+ u = CUU GUU + CU GU + CP GP − − − G− u = CU GU + CP GP + CD GD
For A > 0,
+ =− CUU
+ + where CP+ = 1 − (CUU + CU )
(36)
− − where CP− = 1 − (CU + CD )
hU , 2hUU
+ CU =1+
hU , 2hUU
(37)
− CU = 0.5,
− CD = 0,
Thus, Gu on cell face u in (34) can be expressed using (36) and (37) as: Gu = CUU GUU + CU GU + CP GP + CD GD For A > 0,
CUU = −
hU , 4hUU
CU =
(38)
3 hU , + 4 4hUU
CP =
1 , 4
CD = 0
Similarly, Gd on the cell face d in (34) can be expressed as: Gd = CU GU + CP GP + CD GD + CDD GDD For
A < 0,
CU = 0
CP =
n−1/2 + )
The analytic solutions of (Gp
1 4
n−1/2
3 hD + 4 4hDD
n−1/2 − )
and (Gp
1 3 (Gn−1/2 )+ = − Gn−1 + Gnp , p 2 p 2 Gp
CD =
(39) CDD = −
hD 4hDD
are given as:
(Gn−1/2 )− = p
1 n 1 n+1 G + G 2 p 2 p
(40)
in (35) is obtained using (40) as:
Gn−1/2 = − 41 Gn−1 + Gnp + 41 Gn+1 p p p
(41)
Gs and Gn on the cell face can be expressed similarly. The scheme used here uses a five-point stencil technique. (33) is modified to easily obtain the converged solution. ad au an as GP = Cp (42) GD + GU + GN + GS − S + S Cpx Cpx Cpy Cpy where
S = Cp
1 n−1/2 ad au an as G + (Gd − GD ) + (Gu − GU ) + (Gn − GN ) + (Gs − GS ) . Cpt p Cpx Cpx Cpy Cpy
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
69
5. Computational example A three-dimensional incompressible RANS solver and an unsteady structured grid generation solver, were used in this work to solve three-dimensional problems [13,14]. The Baldwin–Lomax model was used with the wall function to calculate the eddy viscosity in the turbulent flow. In this work, a new computational code for the kinematic boundary condition is developed using a higher order finite analytic scheme. The matrix consisting of the coefficients resulting from a higher order finite analytic scheme was solved using the Generalized Minimal RESidual (GMRES) method [16]. The code is added to the flow solver. Three-dimensional moving free surface turbulent flow around the Wigley hull was computed. The upper boundary moves arbitrarily with the flow, and the grid in the computational domain was generated every time step until the steady state solution was obtained. The explicit schemes are not competitive to solve free surface flow around a ship hull because the time step is chosen sufficiently small (0.0005). Implicit schemes are well known for the property of allowing numerical stability in presence of time steps not restricted by the CFL condition but not free from numerical oscillations. It is usually known that large CFL values leas to instabilities. In this work, all computations were performed using an implicit method without adding any artificial dissipative terms. For unsteady flows, 0.005 time increments were used, which are relatively large. The scheme is able to cope with large CFL numbers. Generally, the skew angle between the velocity components and the faces of the computational cell may lead to increased numerical diffusion. However, the velocity components in the present coordinate system have the same direction as the direction of each curvilinear coordinate line and have physical values. All computations for RANS equations, the pressure Poisson equation and the kinematic boundary condition were performed using a relaxation factor of 0.99, which is more efficient because of its faster convergence rate. With the current methodologies, we hoped to achieve an accurate prediction of the wave patterns around a ship’s hull. The Reynolds number and the Froude number used in the computations and experiment are 3.3 × 106 and 0.316, respectively [10]. No acceleration of the flow was used to reach the final speed. All computations were executed in a grid, sized 127 × 35 × 44. The Wigley hull was used as a ship model. The simulation was performed for the hull advancing in calm water. The example shown here is for the free-surface waves around a Wigley hull. This is an example of an open boundary problem with far field boundaries. The equation for the hull geometry is given by 2 y 2 B 2x 1− z= 1− 2 L D where B is the breath, L is the length, and D id the draft of the ship. The computational domain is one length on the side and one length below the ship hull. The upstream and downstream boundaries were placed at a distance far away from the hull, x/L = −3 and 6, respectively. Waves generated in the interior of the computational domain propagate to the outer boundary and must leave through it. A non-reflecting boundary condition should be applied to pass the wave through the outer boundary. However, this is not a simple problem. The direction of information depends on the sign of A or B in (18). The information comes from inside the computational domain for a positive value of B and outside of the computational domain for a negative value. Zero extrapolation can be used on this boundary for a positive value of B, but for a negative value, the value outside of the computational domain is unknown. Many studies have been done on the appropriate boundary condition to pass the wave through the open boundary. In the present work, a simple Neumann condition was applied at this boundary. The exit plane was set to x = 6. Fortunately, the sign of A in (18) is always positive on the exit plane. Therefore, zero extrapolation can be used on this boundary. The computations were performed for three cases that were identical in all respects except the treatment of the kinematic boundary condition. • Case 1: use a finite analytic scheme based on the 1 − D analytic solution in (29) and allows particles on the upper surface to move only vertically. • Case 2: uses the higher order scheme in (42) and allows particles on the upper surface to move only vertically. • Case 3: uses the higher order scheme in (42) and allows particles on the upper surface to move along the generated curved coordinates lines. For Case 1, only one iteration for each time step was used. Two iterations per time step were used for considering higher order effects. The wave elevation near the hull surface developed quickly by time t = 2 (400 time steps).
70
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
Fig. 4. Comparison of free surface profiles on hull surface, Fn = 0.316.
However, all computations were performed for 1200 time steps (to t = 6). The bow wave height and wave pattern can be used for ranking design variations in ship design. Fig. 4 shows a comparison of the free surface profiles on the hull surface achieved in the three different cases. The two different finite analytic schemes, a lower order and a higher order, were compared. A deviation of wave profiles on the hull surface is observed in the bow region, while better agreement is seen towards the stern. The improvement in the prediction of the free surface profiles on the hull surface was demonstrated by taking into account higher order effects. Figs. 5–7 show wave pattern contours. The two different numerical schemes were also compared. The divergent and the transverse waves are observed. The wave pattern prediction in Fig. 5 is poor especially in the stern region. The higher order scheme was also applied in Figs. 6 and 7. The influence of the higher order effects is found in the far field as well as the near-field waves. The improvement achieved by the present scheme was obtained by using a higher order finite analytic scheme. The particles on the upper surface were allowed to move along the generated curved coordinate lines in Fig. 7. The grid was moved along the hull surface. It was found that there is little difference between the two approaches in (15) and (16). Thus, the present methodologies are suitable for large free surface distortions. (16) was validated by the computations.
Fig. 5. Wave elevation contours (Finite analytic scheme).
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
71
Fig. 6. Wave elevation contours (Higher-order finite analytic scheme, use (15)).
Fig. 7. Wave elevation contours (Higher-order finite analytic scheme, use (16)).
6. Conclusion The kinematic boundary condition was derived by allowing the particles on the upper boundary to move along the generated curved coordinate lines. The new form of the kinematic boundary condition was applied to the computation of free surface flow. This form can be applied to computations around actual ship hulls with a steep grade by allowing the grid to move along the hull surface. A novel approach based on a higher order finite analytic scheme with nonuniform grid spacings was used. In this work, based on the idea of the modified PDEs, a new approach for a hyperbolic PDE was proposed. The results of numerical tests indicated that the method is quite successful. The utility of a higher order finite analytic scheme is demonstrated through the solution of the kinematic boundary condition. It was shown that the methodologies developed in this research contributed to achievement of a sufficiently accurate prediction of the flow solution with the free surface. Appendix A. Supplementary data Supplementary data associated doi:10.1016/j.matcom.2005.11.011.
with
this
article
can
be
found,
in
the
online
version,
at
References [1] M. Beddhu, M.Y. Jiang, L.K. Taylor, D.L. Whitfield, Computation of steady and unsteady flows with a free surface around the Wigley hull, Appl. Math. Comput. 89 (1998) 67–84.
72
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
[2] C.J. Chen, H.C. Chen, Finite analytic method, IIHR Report No. 232-IV. University of Iowa, 1982. [3] C.J. Chen, H.C. Chen, Finite analytic numerical method for unsteady two-dimensional Navier–Stokes equations, J. Comput. Phys. 53 (1984) 209–226. [4] A. Constantin, J. Escher, Symmetry of steady periodic surface water waves with vorticity, J. Fluid Mech. 498 (2004) 171–181. [5] A. Constantin, W. Strauss, Exact steady periodic water waves with vorticity, Comm. Pure Appl. Math. 57 (2004) 481–527. [6] J. Farmer, L. Martinelli, A. Jameson, Fast multigrid method for solving incompressible hydrodynamic problems with free surface, AIAA J. 32 (1994) 1175–1182. [7] J. Farmer, L. Martinelli, G. Cowles, A. Jameson, Fully non-linear CFD techniques for ship performance analysis and design, AIAA (1995) 95–1690. [8] B.R. Hodges, R.L. Street, On simulation of turbulent nonlinear free-surface flows, J. Comput. Phys. 151 (1999) 425–457. [9] ITTC 23rd ITTC Resist. Committee Report, Venuce, Italy, 2002. [10] H. Kajitani, H. Miyata, M. Ikehata, H. Tanaka, H. Adachi, M. Namimatsu, S. Ogiwara, The summary of the cooperative experiment on Wigley parabolic model in Japan, Proeedings of the Second DTRC Workshop on Ship Wave-Resistance Computations, Maryland, 1983. [11] S. Kalliadasis, E.A. Demekhin, C. Ruyer-Quil, M.G. Velarde, Thermocapillary instability and wave formation on a film falling down a uniformly heated plane, J. Fluid Mech. 492 (2003) 303–338. [12] L. Larsson, F. Stern, V. Bertram, Summary and assessment of the Gothburg 2000 workshop, Proceedings of the Gothenburg 2000—A Workshop on Numerical Ship Hydrodynamics, Chalmers University of Technology, Sweden, 2002. [13] S.H. Lee, Three-dimensional uncompressible viscous solutions based on the physical curvilinear coordinate system, PhD Thesis, Mississippi State University, 1997. [14] S.H. Lee, B.K. Soni, The enhancement of an elliptic grid using appropriate control functions, Appl. Math. Comput. 159 (2004) 809–821. [15] S.H. Lee, B.K. Soni, Finite analytic boundary condition for a higher-order finite analytic scheme, Int. J. Numer. Meth. Fluids 50 (3) (2006) 377–396. [16] T. Saad, M.H. Schultz, GMRES: a generalized minimal residual algorithm for solving non-symmetric linear systems, SIAM J. Sci. Stat. Comp. 7 (1986) 856–869. [17] L. Shen, D.K.P. Yue, Large-eddy simulation of free-surface turbulence, J. Fluid Mech. 440 (2001) 75–116. [18] Y. Tahara, F. Stern, An interactive approach for calculating ship boundary layers and wakes for nonzero Froude number, J. Comput. Phys. 98 (1992) 33–53. [19] Y. Tahara, F. Stern, A large-domain approach for calculating ship boundary layers and wakes fields for nonzero Froude number, J. Comput. Phys. 127 (1996) 398–411. [20] C. Truesdell, The physical components of vectors and tensors, J. Appl. Math. Mech. 33 (1953) 345–356. [21] Z.U.A. Warsi, Operations on thephysical components of tensors, ZAMM 76 (1996) 361–362. [22] S. Watanabe, V. Putkaradze, T. Bohr, Integral methods for shallow free-surface flows with separation, J. Fluid Mech. 480 (2003) 233–265.
The derivation and the computation of kinematic boundary condition S.H. Leea,b,∗ , B.K. Sonia,b a
Hydrodynamics Research Department, Hyundai Heavy Industries, Co., Ltd., Ulsan, Korea b Mechanical Engineering, University of Alabama at Birmingham, USA
Received 26 July 2005; received in revised form 27 October 2005; accepted 16 November 2005 Available online 27 December 2005
Abstract The generation of waves by an advancing ship is studied using a moving grid approach. The kinematic boundary condition in the physical curvilinear coordinate system is derived by allowing the particle on the upper boundary to move along the generated curved coordinate lines of the flow field. Thus, the grid is moved along the hull surface. An approach based on the idea of the modified partial differential equation was proposed to discretize a hyperbolic partial differential equation. The computational results obtained using a standard finite analytic scheme were compared with the results obtained using a higher order finite analytic scheme. Convergent steady-state wave profiles for the Wigley hull were compared with experimental measurements. © 2005 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: A moving grid approach; Kinematic boundary condition; Hyperbolic partial differential equation; Higher order finite analytic scheme
1. Introduction In incompressible viscous fluid flows, density is constant and pressure is considered to be a primitive variable. The difficulty in simulating incompressible flow fields lies in solving the coupling of velocity and pressure. Free surface prediction around the ship hull and achieving far-field predictions of sufficient accuracy are still a difficult task [12]. In the proceedings of the 23rd international towing tank conference (ITTC), it was reported that free surface fitting techniques can potentially be quite accurate but are less suitable for large free surface distortions, or for situations in which the grid has to be moved along walls of a complicated shape [9]. These illustrate the basic motivation to resolve the problems in numerical ship hydrodynamics. There are two different methods of treating free surfaces, which depend on the grid type: either fixed or moving. In this work, an approach based on a moving grid was applied to solving the free surface flow. The kinematic boundary condition was then solved with the momentum equations in a coupled manner. A three-dimensional incompressible Reynolds-averaged Navier–Stokes (RANS) solver and an unsteady structured grid generation solver were used in this work [13,14]. A new computational code for the kinematic boundary condition is developed. The code is added to the above flow solver. Three-dimensional moving free surface turbulent flow around the Wigley hull is computed.
∗
Corresponding author. Tel.: +82 52 230 5578; fax: +82 52 230 3410. E-mail address: [email protected] (S.H. Lee)
0378-4754/$32.00 © 2005 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2005.11.011
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
63
The kinematic boundary condition was used to update the unknown location of the free surface. The equation of the free surface F (x, y, z, t) was decomposed into a z-coordinate and G(x, y, t). Particles on the upper boundary of the free surface are assumed to move up and down vertically. This assumption is widely used in wave simulation [8,17,11,22,4]. In the computation of free surface flow around an actual ship hull, the geometry of the hull has a fully complex shape near the fore and stern. The stern region near the upper surface has a steep grade. The grid in this region may have a very large skew, and may make it difficult to obtain a reasonable solution with this assumption. To avoid this problem, this work newly derives the kinematic boundary condition in the physical curvilinear coordinate system by allowing the particle on the upper boundary to move along the generated curved coordinate lines in the flow field. In the physical curvilinear coordinate system, each component of the velocity vector has the same direction as that of each curvilinear coordinate line and has a physical value. The physical curvilinear component form of the velocity was first introduced by Truesdell [20]. Lee [13] and Warsi [21] derived formulas needed in order to transform differential equations into a physical curvilinear coordinate system. In this work, Lee [13] and Warsi [21]’s method is used. Another important issue is the discretization of the kinematic boundary condition. This condition is a hyperbolic partial differential equation (PDE). Several approaches have been tried to obtain a reasonable discretization. Tahara and Stern [18] used a backward difference for the x-derivative term and a central difference scheme for the y-derivative term. Later, Tahara and Stern [19] used a Beam and Warming linear-multi-step method. They added implicit and explicit fourth-order artificial dissipations to a difference scheme. Beddhu et al. [1] used an explicit scheme to solve this condition without adding an artificial dissipative term. Farmer et al. [6] used a local time stepping scheme and added background dissipation to prevent decoupling of the solution. Later, Farmer et al. [7] used symmetric-limited positive and upstream-limited positive schemes to avoid overturning of the waves. In this work, an approach based on analytic methods is used to discretize the kinematic boundary condition, which requires no artificial dissipative terms. Recently, Lee and Soni [15] derived higher-order finite analytic schemes based on one-dimensional local analytic solutions to solve three-dimensional PDEs. In the present research, the same higher order schemes were derived and applied to two-dimensional hyperbolic PDEs using the idea of the modified PDE. 2. Physical curvilinear coordinate system In the following analysis, xk are Cartesian coordinates, ξ i are curvilinear coordinates and ξ (i) are physical curvilinear coordinates. The relationship between the coordinates ξ i and ξ (i) is defined as follows: √ ξ (i) = ξ (i) (ξ i ) where ξ (i) = gii ξ i , ξ 1 = ξ 2 = ξ 3 = 1 (1) The coordinates ξ (i) resemble the coordinate stretching in each direction of ξ i . In order to transform the coordinates, the vector dr can be written as: dr = ai dξ i = a(i) dξ (i)
where ai =
∂r ∂r and a(i) = (i) ∂ξ i ∂ξ
(2)
ai are covariant base vectors and a(i) are covariant base vectors in the physical curvilinear coordinate system. The relationship between the two covariant base vectors is obtained from (2) as: 1 a(i) = √ ai gii
(3)
The physical curvilinear components of a velocity vector, u(i) , are defined as the magnitude of the ith component projected onto the ith physical curvilinear coordinate direction: u(i) = u · a(i)
(4)
This form was introduced to give physical meaning to the contravariant component of a velocity lost by full transformation. The relationship between the velocity components in the two different coordinate systems is defined as: u(i) = u(j) a(j) · ei = u(j)
∂xi ∂xi 1 = √ u(j) j ∂ξ (j) gjj ∂ξ
(5)
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
64
u(j) = u(i)ei · a(j) = u(i)
∂ξ (j) ∂ξ j √ = gjj u(i) ∂xi ∂xi
(6)
The velocity components in u(i) in Eq. (5) represent (u, v, w). The repeated indices above imply a sum. Using the general transformation laws, for a scalar φ, the gradient is written as: ∇φ =
∂φ (i) ∂φ a = g(ik) (k) a(i) ∂ξ (i) ∂ξ
where a(i) = g(ik) a(k)
and a(i) =
(7)
√ gii ai
ai are contravariant base vectors and a(i) are contravariant base vectors in the physical curvilinear coordinate system. The relationship between the two contravariant metric tensors in (7) is written as: g(ij) =
√ √ gii gjj gij
where gij =
∂ξ i ∂ξ j ∂xk ∂xk
(8)
Also, the relationship between two covariant metric tensors is written as: g(ij) = √
1 √
gii gjj
gij
where gij =
∂xk ∂xk ∂ξ i ∂ξ j
(9)
g(ij) and g(ij) above are the covariant and contravariant metric tensors in the physical curvilinear coordinate system, respectively. The dot product of the two vectors u and ∇φ is expressed as: u · ∇φ = u(i) a(i) ·
∂φ (j) ∂φ a = u(i) (i) ∂ξ (j) ∂ξ
(10)
3. Kinematic boundary condition The concept of the material surface has been used to derive the kinematic boundary condition, which expresses the fact that the free surface moves with the fluid and deforms in shape as the flow progresses. There is no mass transport through the surface, so its boundary is always composed of the same fluid particles. If the equation of the surface is represented by F = 0, then the equation of the surface at t + t is expressed as: F (r + δr, t + t) = 0
(11)
Using a Taylor’ series expansion, the kinematic boundary condition can be derived by taking the first order components in the vector form: ∂F + u · ∇F = 0 (12) ∂t It is assumed here that flow particles move up and down vertically. If the equation of the surface F (x, y, z, t) is decomposed into the z-coordinate and G(x, y, t), (12) can be written in the Cartesian coordinates system by substituting the z-coordinate and G for the surface, F as follows: ∂G ∂G ∂G +u +v =w ∂t ∂x ∂y
where F (x, y, z, t) = z − G(x, y, t)
(13)
(12) can also be expressed in the physical curvilinear coordinate system. The vector form of (12) is transformed by using (10) as follows: ∂F ∂F + u(i) (i) = 0 ∂t ∂ξ
(14)
If we assume that the flow particles move only vertically, (14) can be rewritten as follows: ∂G ∂G ∂G + u(1) (1) + u(2) (2) = w ∂t ∂ξ ∂ξ
where w = u(i)
∂z ∂ξ (i)
(15)
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
65
Fig. 1. The particles on the upper boundary moving up and down vertically.
Fig. 2. The particles on the upper boundary moving along the generated curved coordinate lines.
Here, the surface equation F was decomposed into its z-coordinate and G(ξ (1) , ξ (2) , t). (15) can be used where the vertical movement is dominant [5]. Fig. 1 shows the particles on the upper boundary moving up and down vertically. When flow particles move along the generated curved coordinate lines ξ (3) , (14) can be rewritten as follows: ∂G ∂G ∂G ∂H ∂H + u(1) (1) + u(2) (2) = u(1) (1) + u(2) (2) + u(3) ∂t ∂ξ ∂ξ ∂ξ ∂ξ
(16)
The surface equation F was decomposed into H(ξ (3) ) and G(ξ (1) , ξ (2) , t). H represents the line in the ξ (3) -coordinate direction. Fig. 2 shows the particles on the upper boundary moving along the generated curved coordinate lines. 4. Discretization methods 4.1. Finite analytic scheme The kinematic boundary condition is coupled with momentum equations. This boundary condition is required to regenerate the grid in the flow fields at every time step. Thus, it is necessary to solve this condition accurately. The scheme used here can be applied to the discretization of hyperbolic PDEs. (16) are rearranged to make the standard form of the finite analytic method developed by Chen and Chen [3]. 0=a
∂G ∂G ∂G + b (2) + c + S ∂ξ (1) ∂ξ ∂t
(17)
where S = −u(1)
∂H ∂H − u(2) (2) − u(3) , ∂ξ (1) ∂ξ
a = u(1) , b = u(2) , c = 1
The second derivative terms are added to the right hand side of (17). To minimize the effects of the second derivative terms, each second derivative term can be multiplied by a very small value. Thus modified, the equation for the kinematic boundary condition can be written as follows: ∂2 G ∂ξ
(1) 2
+
∂2 G ∂ξ
(2) 2
+
∂2 G ∂G ∂G ∂G = A (1) + B (2) + C +S ∂ξ ∂ξ ∂t ∂t 2
(18)
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
66
where
1 (1) ∂H (2) ∂H (3) S= −u , +u −u α ∂ξ (1) ∂ξ (2)
A=
u(1) , α
B=
u(2) , α
C=
1 α
(18) is the standard form of the finite analytic method that solves the elliptic PDE. A and B are assumed to be constant over the local element. A six-point finite analytic scheme based on a one-dimensional analytic solution is derived at the vertex of a cell. (18) is rearranged by decomposing into one-dimensional equations as shown below: ∂2 G 2 ∂ξ (1)
∂2 G 2 ∂ξ (2)
−A
∂G = S1 ∂ξ (1)
(19)
−B
∂G = S2 ∂ξ (2)
(20)
∂G ∂2 G −C = S3 2 ∂t ∂t
(21)
S1 + S2 + S3 = S
(22)
where S1 = − S3 = −
∂2 G ∂ξ (2)
2
∂2 G ∂ξ
(1) 2
+B
∂G ∂2 G ∂G − +C + S, 2 (2) ∂ξ ∂t ∂t
+A
∂G ∂2 G ∂G − + B (2) + S 2 (1) (2) ∂ξ ∂ξ ∂ξ
S2 = −
∂2 G ∂ξ (1)
2
+A
∂G ∂2 G ∂G − +C + S, 2 (1) ∂ξ ∂t ∂t
The source terms S1 , S2 and S3 are assumed to be constant in a local cell. The analytic solution for one-dimensional equation [2] is defined in the present coordinate system using three nodal points. The solution for (19)–(21) is expressed as: (1)
− 1) + b1 −
S1 A
(23)
(2)
− 1) + b2 −
S2 B
(24)
G(ξ (1) ) = a1 (eAξ G(ξ (2) ) = a2 (eBξ
G(t) = a3 (eCt − 1) + b3 −
S3 C
(25)
The coefficients a1 , b1 , a2 , b2 , a3 and b3 are determined by applying the boundary conditions at ξ (1) = hD and −hU , ξ (2) = hN and −hS , respectively. The subscripts U and D denote upstream and downstream, respectively. Fig. 3 shows local elements. The finite analytic solution for the one-dimensional modified PDE can be obtained by taking sufficiently large values of A, B and C as: GP = aD GD + aU GU − CPX S1 If A > 0,
aD = 0, aU = 1, CPX =
(26) hU A
else
aD = 1, aU = 0, CPX =
−hD A
GP = aN GN + aS GS − CPY S2 If B > 0,
aN = 0, aS = 1, CPY =
GP = aT Gn+1 + aB GnP − CPT S3 P aT = 0, aB = 1, CPT =
t C
(27) hS B
else
aN = 1, aS = 0, CPY =
−hN B (28)
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
67
Fig. 3. Local elements.
A six-point finite analytic scheme is derived by substituting the equations for S1 , S2 and S3 from (26)–(28) into (22) as: 1 aD aU aN aS n GP = CP G + GD + GU + GN + GS − S (29) CPT P CPX CPX CPY CPY where 1 1 1 1 = + + CP CPX CPY CPT 4.2. Higher-order finite analytic scheme The one-dimensional equation in (19)–(21) can be solved analytically in a local cell by using one center point and two end cell faces as: Gp = ad Gd + au Gu − Cpx S1
(30)
Gp = an Gn + as Gs − Cpy S2
(31)
Gp = at Gn+1/2 + ab Gn−1/2 − Cpt S3 p p
(32)
ad = at =
1 − e−Ahu , eAhd − e−Ahu 1 − e−Ct/2 , − e−Ct/2
eCt/2
au =
eAhd − 1 , eAhd − e−Ahu
aw =
an =
eCt/2 − 1 , − e−Ct/2
eCt/2
1 − e−Bhs , eBhn − e−Bhs Cpx =
as =
eBhn − 1 , eBhn − e−Bhs
1 hd (e−Ahu − 1) + hu (eAhd − 1) , A eAhd − e−Ahu
1 hn (e−Bhs − 1) + hs (eBhn − 1) t (e−Ct/2 − 1) + (eCt/2 − 1) , C = pt B eBhn − e−Bhs 2C eCt/2 − e−Ct/2 Here, the subscripts u, d, n and s on G represent the values on the cell faces. The value at point p is obtained from the values at the cell faces d, u, n and s. Fig. 3 shows the grid size of a local cell. For sufficiently large values of A, B and C, the analytic coefficients in (30)–(32) can be approximated as: Cpy =
hu −hd else ad = 1, au = 0, Cpx = , 2A 2A hs −hn t = else an = 1, as = 0, Cpy = , at = 0, ab = 1, Cpt = 2B 2B 2C
If A > 0,
ad = 0, au = 1, Cpx =
If B > 0,
an = 0, as = 1, Cpy
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
68
A formula for an unsteady two-dimensional equation is obtained as: ab n−1/2 ad au an as Gp = Cp Gp + Gd + Gu + Gn + Gs − S Cpt Cpx Cpx Cpy Cpy
(33)
Here, the value on cell face u in (33) can be solved analytically for a local cell by applying the boundary condition at points U and D, or UU and P. However, the value G− u on cell face u obtained by applying the boundary condition at points U and D is different from the value G+ u obtained at points UU and P. Therefore, in this work, G on the cell face in (33) will be defined as [15]: Gd =
− G+ d + Gd , 2
Gn−1/2 = p
n−1/2 + )
(Gp
Gu =
− G+ u + Gu , 2
Gn =
− G+ n + Gn , 2
Gs =
− G+ s + Gs 2
n−1/2 − )
+ (Gp 2
(34) (35)
− The analytic solutions for G+ u and Gu in (34) are obtained by applying the boundary condition at GUU and GU , GU and GD , respectively: + + + G+ u = CUU GUU + CU GU + CP GP − − − G− u = CU GU + CP GP + CD GD
For A > 0,
+ =− CUU
+ + where CP+ = 1 − (CUU + CU )
(36)
− − where CP− = 1 − (CU + CD )
hU , 2hUU
+ CU =1+
hU , 2hUU
(37)
− CU = 0.5,
− CD = 0,
Thus, Gu on cell face u in (34) can be expressed using (36) and (37) as: Gu = CUU GUU + CU GU + CP GP + CD GD For A > 0,
CUU = −
hU , 4hUU
CU =
(38)
3 hU , + 4 4hUU
CP =
1 , 4
CD = 0
Similarly, Gd on the cell face d in (34) can be expressed as: Gd = CU GU + CP GP + CD GD + CDD GDD For
A < 0,
CU = 0
CP =
n−1/2 + )
The analytic solutions of (Gp
1 4
n−1/2
3 hD + 4 4hDD
n−1/2 − )
and (Gp
1 3 (Gn−1/2 )+ = − Gn−1 + Gnp , p 2 p 2 Gp
CD =
(39) CDD = −
hD 4hDD
are given as:
(Gn−1/2 )− = p
1 n 1 n+1 G + G 2 p 2 p
(40)
in (35) is obtained using (40) as:
Gn−1/2 = − 41 Gn−1 + Gnp + 41 Gn+1 p p p
(41)
Gs and Gn on the cell face can be expressed similarly. The scheme used here uses a five-point stencil technique. (33) is modified to easily obtain the converged solution. ad au an as GP = Cp (42) GD + GU + GN + GS − S + S Cpx Cpx Cpy Cpy where
S = Cp
1 n−1/2 ad au an as G + (Gd − GD ) + (Gu − GU ) + (Gn − GN ) + (Gs − GS ) . Cpt p Cpx Cpx Cpy Cpy
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
69
5. Computational example A three-dimensional incompressible RANS solver and an unsteady structured grid generation solver, were used in this work to solve three-dimensional problems [13,14]. The Baldwin–Lomax model was used with the wall function to calculate the eddy viscosity in the turbulent flow. In this work, a new computational code for the kinematic boundary condition is developed using a higher order finite analytic scheme. The matrix consisting of the coefficients resulting from a higher order finite analytic scheme was solved using the Generalized Minimal RESidual (GMRES) method [16]. The code is added to the flow solver. Three-dimensional moving free surface turbulent flow around the Wigley hull was computed. The upper boundary moves arbitrarily with the flow, and the grid in the computational domain was generated every time step until the steady state solution was obtained. The explicit schemes are not competitive to solve free surface flow around a ship hull because the time step is chosen sufficiently small (0.0005). Implicit schemes are well known for the property of allowing numerical stability in presence of time steps not restricted by the CFL condition but not free from numerical oscillations. It is usually known that large CFL values leas to instabilities. In this work, all computations were performed using an implicit method without adding any artificial dissipative terms. For unsteady flows, 0.005 time increments were used, which are relatively large. The scheme is able to cope with large CFL numbers. Generally, the skew angle between the velocity components and the faces of the computational cell may lead to increased numerical diffusion. However, the velocity components in the present coordinate system have the same direction as the direction of each curvilinear coordinate line and have physical values. All computations for RANS equations, the pressure Poisson equation and the kinematic boundary condition were performed using a relaxation factor of 0.99, which is more efficient because of its faster convergence rate. With the current methodologies, we hoped to achieve an accurate prediction of the wave patterns around a ship’s hull. The Reynolds number and the Froude number used in the computations and experiment are 3.3 × 106 and 0.316, respectively [10]. No acceleration of the flow was used to reach the final speed. All computations were executed in a grid, sized 127 × 35 × 44. The Wigley hull was used as a ship model. The simulation was performed for the hull advancing in calm water. The example shown here is for the free-surface waves around a Wigley hull. This is an example of an open boundary problem with far field boundaries. The equation for the hull geometry is given by 2 y 2 B 2x 1− z= 1− 2 L D where B is the breath, L is the length, and D id the draft of the ship. The computational domain is one length on the side and one length below the ship hull. The upstream and downstream boundaries were placed at a distance far away from the hull, x/L = −3 and 6, respectively. Waves generated in the interior of the computational domain propagate to the outer boundary and must leave through it. A non-reflecting boundary condition should be applied to pass the wave through the outer boundary. However, this is not a simple problem. The direction of information depends on the sign of A or B in (18). The information comes from inside the computational domain for a positive value of B and outside of the computational domain for a negative value. Zero extrapolation can be used on this boundary for a positive value of B, but for a negative value, the value outside of the computational domain is unknown. Many studies have been done on the appropriate boundary condition to pass the wave through the open boundary. In the present work, a simple Neumann condition was applied at this boundary. The exit plane was set to x = 6. Fortunately, the sign of A in (18) is always positive on the exit plane. Therefore, zero extrapolation can be used on this boundary. The computations were performed for three cases that were identical in all respects except the treatment of the kinematic boundary condition. • Case 1: use a finite analytic scheme based on the 1 − D analytic solution in (29) and allows particles on the upper surface to move only vertically. • Case 2: uses the higher order scheme in (42) and allows particles on the upper surface to move only vertically. • Case 3: uses the higher order scheme in (42) and allows particles on the upper surface to move along the generated curved coordinates lines. For Case 1, only one iteration for each time step was used. Two iterations per time step were used for considering higher order effects. The wave elevation near the hull surface developed quickly by time t = 2 (400 time steps).
70
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
Fig. 4. Comparison of free surface profiles on hull surface, Fn = 0.316.
However, all computations were performed for 1200 time steps (to t = 6). The bow wave height and wave pattern can be used for ranking design variations in ship design. Fig. 4 shows a comparison of the free surface profiles on the hull surface achieved in the three different cases. The two different finite analytic schemes, a lower order and a higher order, were compared. A deviation of wave profiles on the hull surface is observed in the bow region, while better agreement is seen towards the stern. The improvement in the prediction of the free surface profiles on the hull surface was demonstrated by taking into account higher order effects. Figs. 5–7 show wave pattern contours. The two different numerical schemes were also compared. The divergent and the transverse waves are observed. The wave pattern prediction in Fig. 5 is poor especially in the stern region. The higher order scheme was also applied in Figs. 6 and 7. The influence of the higher order effects is found in the far field as well as the near-field waves. The improvement achieved by the present scheme was obtained by using a higher order finite analytic scheme. The particles on the upper surface were allowed to move along the generated curved coordinate lines in Fig. 7. The grid was moved along the hull surface. It was found that there is little difference between the two approaches in (15) and (16). Thus, the present methodologies are suitable for large free surface distortions. (16) was validated by the computations.
Fig. 5. Wave elevation contours (Finite analytic scheme).
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
71
Fig. 6. Wave elevation contours (Higher-order finite analytic scheme, use (15)).
Fig. 7. Wave elevation contours (Higher-order finite analytic scheme, use (16)).
6. Conclusion The kinematic boundary condition was derived by allowing the particles on the upper boundary to move along the generated curved coordinate lines. The new form of the kinematic boundary condition was applied to the computation of free surface flow. This form can be applied to computations around actual ship hulls with a steep grade by allowing the grid to move along the hull surface. A novel approach based on a higher order finite analytic scheme with nonuniform grid spacings was used. In this work, based on the idea of the modified PDEs, a new approach for a hyperbolic PDE was proposed. The results of numerical tests indicated that the method is quite successful. The utility of a higher order finite analytic scheme is demonstrated through the solution of the kinematic boundary condition. It was shown that the methodologies developed in this research contributed to achievement of a sufficiently accurate prediction of the flow solution with the free surface. Appendix A. Supplementary data Supplementary data associated doi:10.1016/j.matcom.2005.11.011.
with
this
article
can
be
found,
in
the
online
version,
at
References [1] M. Beddhu, M.Y. Jiang, L.K. Taylor, D.L. Whitfield, Computation of steady and unsteady flows with a free surface around the Wigley hull, Appl. Math. Comput. 89 (1998) 67–84.
72
S.H. Lee, B.K. Soni / Mathematics and Computers in Simulation 71 (2006) 62–72
[2] C.J. Chen, H.C. Chen, Finite analytic method, IIHR Report No. 232-IV. University of Iowa, 1982. [3] C.J. Chen, H.C. Chen, Finite analytic numerical method for unsteady two-dimensional Navier–Stokes equations, J. Comput. Phys. 53 (1984) 209–226. [4] A. Constantin, J. Escher, Symmetry of steady periodic surface water waves with vorticity, J. Fluid Mech. 498 (2004) 171–181. [5] A. Constantin, W. Strauss, Exact steady periodic water waves with vorticity, Comm. Pure Appl. Math. 57 (2004) 481–527. [6] J. Farmer, L. Martinelli, A. Jameson, Fast multigrid method for solving incompressible hydrodynamic problems with free surface, AIAA J. 32 (1994) 1175–1182. [7] J. Farmer, L. Martinelli, G. Cowles, A. Jameson, Fully non-linear CFD techniques for ship performance analysis and design, AIAA (1995) 95–1690. [8] B.R. Hodges, R.L. Street, On simulation of turbulent nonlinear free-surface flows, J. Comput. Phys. 151 (1999) 425–457. [9] ITTC 23rd ITTC Resist. Committee Report, Venuce, Italy, 2002. [10] H. Kajitani, H. Miyata, M. Ikehata, H. Tanaka, H. Adachi, M. Namimatsu, S. Ogiwara, The summary of the cooperative experiment on Wigley parabolic model in Japan, Proeedings of the Second DTRC Workshop on Ship Wave-Resistance Computations, Maryland, 1983. [11] S. Kalliadasis, E.A. Demekhin, C. Ruyer-Quil, M.G. Velarde, Thermocapillary instability and wave formation on a film falling down a uniformly heated plane, J. Fluid Mech. 492 (2003) 303–338. [12] L. Larsson, F. Stern, V. Bertram, Summary and assessment of the Gothburg 2000 workshop, Proceedings of the Gothenburg 2000—A Workshop on Numerical Ship Hydrodynamics, Chalmers University of Technology, Sweden, 2002. [13] S.H. Lee, Three-dimensional uncompressible viscous solutions based on the physical curvilinear coordinate system, PhD Thesis, Mississippi State University, 1997. [14] S.H. Lee, B.K. Soni, The enhancement of an elliptic grid using appropriate control functions, Appl. Math. Comput. 159 (2004) 809–821. [15] S.H. Lee, B.K. Soni, Finite analytic boundary condition for a higher-order finite analytic scheme, Int. J. Numer. Meth. Fluids 50 (3) (2006) 377–396. [16] T. Saad, M.H. Schultz, GMRES: a generalized minimal residual algorithm for solving non-symmetric linear systems, SIAM J. Sci. Stat. Comp. 7 (1986) 856–869. [17] L. Shen, D.K.P. Yue, Large-eddy simulation of free-surface turbulence, J. Fluid Mech. 440 (2001) 75–116. [18] Y. Tahara, F. Stern, An interactive approach for calculating ship boundary layers and wakes for nonzero Froude number, J. Comput. Phys. 98 (1992) 33–53. [19] Y. Tahara, F. Stern, A large-domain approach for calculating ship boundary layers and wakes fields for nonzero Froude number, J. Comput. Phys. 127 (1996) 398–411. [20] C. Truesdell, The physical components of vectors and tensors, J. Appl. Math. Mech. 33 (1953) 345–356. [21] Z.U.A. Warsi, Operations on thephysical components of tensors, ZAMM 76 (1996) 361–362. [22] S. Watanabe, V. Putkaradze, T. Bohr, Integral methods for shallow free-surface flows with separation, J. Fluid Mech. 480 (2003) 233–265.